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Summary 


We  performed  research  in  the  field  of  electromagnetic  wave  chaos  and  High 
Power  Microwave  (HPM)  effects  in  electronic  systems.  Our  emphasis  was  on  issues 
likely  to  be  most  relevant  to  the  coupling  of  electromagnetic  radiation  into  systems  and 
its  effects  on  systems,  as  in  electromagnetic  attacks  of  military  electronics.. 

Our  program  addressed  the  following  issues  of  relevance  to  the  understanding  of 
HPM  effects  on  electronic  systems: 

•  The  extension  of  the  previously  developed  Random  Coupling  Model  statistical 
description  of  wave  coupling  into  enclosures  to  describe  a)  the  coupling  through 
apertures,  b)  the  coupling  to  mixed  systems  for  which  only  part  of  the  ray  phase  space  is 
chaotic,  and  c)  the  coupling  to  systems  of  systems  in  which  the  components  have  varying 
degrees  of  isolation,  for  example  chains  of  cavities. 

•  The  further  development  of  time  domain  models  for  the  response  of  systems  excited  by 
pulses  of  wave  energy,  including  nonlinear  effects 


•  The  Extension  of  HPM  upset  experimentation  and  modeling  to  complex  networks  of 
interconnected  circuits .  Our  aim  is  to  develop  a  generalized  approach  for  modeling  HPM 
effects  in  electronic  circuits,  systems  and  infrastructures. 

•  Examination  of  fading,  power-delay  profiles,  and  small-signal  discrimination  in 
reverberant  electromagnetic  environments  containing  short-ray  trajectories. 


General  theory  of  HPM  coupling  to  circuits  in  enclosures  through  apertures 

The  RCM  is  based  on  the  use  of  Random  Matrix  Theory  (RMT)  that  has  found 
wide  application  in  mesoscopic  and  nuclear  physics.  The  model  incorporates  system- 
specific  information  about  the  near  field  behavior  of  the  ports  of  the  enclosure,  the 
volume  of  the  enclosure,  and  a  measure  of  the  enclosure’s  average  quality  factor,  and 
predicts  the  statistical  behavior  of  the  enclosures  scattering  parameters  (actually  the 
impedance  matrix  from  which  the  scattering  matrix  can  be  obtained).  This  model  was 
introduced  in  Refs.  1-3  for  a  specific  geometry,  namely  planar  cavities  excited  by  point 
like  sources.  In  Ref  6.  we  developed  a  formulation  of  the  model  for  the  case  of  three- 
dimensional  cavities,  with  arbitrary  polarization  of  the  fields  inside,  with  ports  that  are 
large  compared  with  a  wavelength  and  with  conductors  inside  that  act  as  antennas.  The 
use  of  this  model  allows  one  to  study  the  coupling  process  starting  from  an  incident  plane 
wave  that  penetrates  through  an  aperture,  distributing  its  energy  throughout  an  enclosure 
and  inducing  voltages  on  conductors  in  the  enclosure.  This  is  done  by  breaking  the 
problem  into  pieces.  The  pieces  are  then  combined  in  the  RCM  to  give  a  statistical 
description  of  the  coupling  processes. 
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Figure  1:  Schematic  of  an  enclosure  with  an  aperture  illuminated  from  the  outside  and  an  antenna 
(representing  a  circuit  element)  inside  that  is  connected  to  a  load  ZL.  Coupling  from  the  incident  wave  to 
the  antenna  is  described  by  a  system  of  linear  matrix  equations.  Statistics  (and  the  RCM)  enter  through  the 
'/-matrix  defined  in  the  text. 


Figure  1  illustrates  how  this  is  done.  A  wave  is  incident  on  the  enclosure  illuminating  the 
aperture  and  is  partially  reflected.  The  field  on  the  aperture  is  represented  in  a  basis  of 
modes  with  amplitudes  Vs .  The  properties  of  the  aperture  are  described  by  the  aperture 

radiation  admittance  Ymd  =  /  Im^K""']  +  G""' ,  which  gives  the  amplitudes  Is  =  YrsfVs,  of  the 

magnetic  field  basis  functions  when  the  enclosure  is  absent  and  the  aperture  radiates  into 
free  space.  Similarly  the  currents  in  the  antenna  can  be  represented  in  a  basis  of  modes 
with  amplitudes  I .  The  properties  of  the  antenna  are  described  by  its  radiation 


impedance  Zrarf  =  ;Im[z™'i]  +  ^""/ giving  the  voltage  on  the  antenna  when  radiating  into 
free  space.  When  the  antenna  (attached  to  an  effective  load  ZL )  is  placed  in  the 
enclosure  it  becomes  coupled  to  the  fields  in  the  aperture.  The  precise  characterization  of 
this  coupling  is  a  complicated  computational  electromagnetics  problem.  However,  if  the 
enclosure  has  chaotic  trajectories,  the  relationship  between  antenna  and  aperture  voltages 
and  currents  is  modeled  by  the  matrix  T  =  j  I m |  U ]  +  [V]12  •  <g  •  f V]12  •  Here  the  matrices  U 


and  V  are  given  in  terms  of  the  radiation  admittance  of  the  aperture  and  the  radiation 


impedance  of  the  antenna, 
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All  the  statistics,  and  coupling  between  the  aperture  and  antenna,  are  imbedded  in  the 
random  matrix  E,  whose  properties  are  given  in  terms  of  the  eigenfunctions  and 


eigenvalues  of  a  Guassian  Random  Matrix  and  depend  only  on  the  average  internal  losses 
in  the  cavity  and  the  average  density  of  modes.  The  system  specific  information  is 
provided  by  the  marices  Zmd  and  Ymd  that  can  be  computed  or  measured  in  isolation 

from  the  rest  of  the  system.  This  formulation  provides  a  compact  approach  to  analyzing 
coupling  for  specific  configurations  based  on  the  properties  of  the  components  of  the 
configuration.  It  has  no  adjustable  parameters.  We  are  currently  in  the  process  of  testing 
this  formulation  in  the  laboratory. 

Distribution  of  RF  power  in  chains  of  coupled  cavities 

Our  previous  studies  focused  on  the  distribution  of  RF  energy  throughout  a  single 
enclosed  space.  However,  a  problem  that  is  often  of  interest  is  the  distribution  of  energy 

through  a  chain  of  enclosed  spaces  that 
are  coupled  together.  These  may  be 
either  compartments  in  a  vehicle  or 
aircraft,  or  boxes  enclosed  within  other 
boxes.  We  have  extended  our  RCM  to 
the  study  of  the  problem  of  the  statistics 
of  fields  in  a  chain  of  coupled  cavities 
[13].  In  the  general  case,  the  coupled 
power  is  given  in  terms  of  a  product  of 
random  matrices.  The  distribution  of 
coupled  power  then  has  the  same 
statistics  as  the  finite  time  Lyapunov 
numbers  for  a  discrete  dynamical 
system.  This  distribution  approaches 
log-normal  as  the  number  of  cavities 
(and  matrices)  is  increased.  However, 
there  are  significant  large  deviation  departures  from  log-normal  which  describe  the  rare 
cases  of  large  coupled  power  and  are  thus  of  interest.  We  have  derived  general 
expressions  for  the  probability  distribution  of  the  power  coupled  into  the  Nth  cavity  of  a 
chain  in  the  case  in  which  losses  in  the  cavities  is  the  dominant  absorption  mechanism. 


Figure  2  Distributions  of  power  ratios  for  chains  of 
up  to  seven  cavities  with  moderate  losses  for  a  large 
number  of  cavities  the  distribution  becomes  log- 
normal,  from  [23] 


These  expressions,  which  are  in  the  form  of  products  of  system  specific  factors 
describing  the  coupling  ports  and  losses,  and  universal  statistical  factors,  have  been 
compared  with  Monte  Carlo  evaluations  of  the  coupling  using  the  RCM. 

Given  that  the  power  coupled  through  a  chain  of  cavities  is  a  product  of  factors,  it 
is  natural  that  the  probability  distribution  for  coupled  power  tends  to  log-normal  as  the 
number  of  cavities  increases.  This  is  illustrated  in  Fig.  2  where  the  power  transmission 
factor  ( R[N)  =  power  to  cavity  N  /  power  to  cavity  1)  for  chains  of  up  to  seven  cavities  are 
plotted  for  a  case  in  which  the  cavities  have  moderate  losses  (Q  -  width  equals  mode 
spacing).  The  present  model  assumes  that  successive  cavities  are  coupled  by  ports  or 
transmission  lines  that  support  a  single  mode.  In  the  coming  year  we  intend  to  extend 
this  to  coupling  through  over-moded  or  resonant  apertures. 

Fading  in  high  and  low  loss  environments 

Fading  is  the  observation  of  variations  in  signal  strength  measured  at  a  receiver 
due  to  time-dependent  variations  in  the  propagation  of  waves  from  the  source,  or  due  to 
multi-path  scattering  and  interference.  It  is  well  known  that  the  quantitative  statistical 
theory  of  wave  chaos  -  random  matrix  theory  (RMT)  -  can  be  successfully  applied  to 
predict  statistical  properties  of  many  quantities,  such  as  the  scattering  matrix,  of  a  wave 
chaotic  system.  We  started  from  the  statistical  model  of  the  scattering  matrix  [14]  to 
establish  a  general  fading  model  valid  under  all  conditions.  The  model  provides  a  first 
principles  understanding  of  the  most  common  statistical  model  used  in  the 
communications  field,  namely  Rayleigh  fading,  and  shows  that  the  statistical  properties 
are  governed  by  a  single  quantity  related  to  the  loss  or  de-phasing  parameter  of  RMT. 
We  also  combined  the  RMT  fading  model  with  our  random  coupling  model  (RCM)  that 
takes  into  account  system-specific  features  such  as  direct  and  short  orbits  [13,  25-27],  to 
build  a  more  general  fading  model  that  includes  Rician  fading. 


Figure  3:  The  probability  density  functions  of  the  fading  amplitude  I  srmt, 21  for  the  time  reversal 

invariant  case  (left)  and  srmt, 21  for  the  broken  time  reversal  invariant  case  (right).  Solid  curves  show  the 
numerical  results  from  the  RMT  model  with  different  loss  parameters  ( O  =  0,  0.1,  1,  10).  For  the  higher 
loss  cases  ( a  =  1  and  a  =  10),  the  corresponding  Rayleigh  distributions  are  shown  as  dashed  curves. 

In  the  high  loss-parameter  limit,  our  model  agrees  with  the  Rayleigh/Rice  models, 
however  it  shows  significant  deviations  from  the  Rayleigh/Rice  distribution  in  the  limit 
of  low  loss  [17].  We  have  performed  experiments  with  two  ray-chaotic  microwave 
cavities  [15-17]  to  test  the  RMT/RCM  fading  model  over  a  wide  range  of  loss  parameter 
values,  and  our  results  are  in  excellent  agreement  with  our  theory. 


Statistics  of  field  distribution  in  cavities  with  both  regular  and  chaotic  orbits 

This  work  was  included  in  the  Ph.D.  thesis  of  Ming-jer  Lee,  who  graduates  8/13. 
The  RCM  is  based  on  idea  that  the  ray  trajectories  describing  the  propagation  of  radiation 
throughout  an  enclosure  are  chaotic.  However,  for  some  enclosures  the  walls  are  flat  and 
facing  each  other  in  such  a  way  that  some  ray  trajectories  are  not  chaotic  while  others  are. 
We  have  extended  our  model  to  an  example  of  such  a  cavity.  We  have  found  that  the 
response  of  the  cavity  to  excitations  in  this  case  can  be  characterized  as  a  sum  of  the 
response  of  modes  whose  trajectories  are  chaotic  and  modes  whose  trajectories  are  not 
chaotic.  The  chaotic  portion  of  the  response  has  the  same  universal  features  that  the 
RCM  predicts.  The  nonchaotic  response  of  the  system  has  different  features  that  can  still 
be  characterized  generally,  that  is,  without  detailed  solution  of  the  field  equations  in  the 
cavity. 

An  example  of  the  effect  that  we  consider  is  illustrated  in  Fig.  4a  and  Fig.  4b. 
Figure  4a  shows  two  eigenmodes  (numbers  5005  and  5006)  of  a  two-dimensional 
mushroom  shaped  cavity  that  have  nearly  the  same  eigenvalue.  The  one  on  the  left 
extends  throughout  the  cavity  and  corresponds  to  rays  that  ergodically  fill  the  cavity, 
while  the  one  on  the  right  extends  mainly  to  the  semicircular  region  and  corresponds  to 
rays  bounce  along  the  vertical,  horizontal,  and  outer  circular  boundary.  The  histogram  for 
impedances  for  a  port  that  is  placed  in  the  upper  right  hand  side  (which  will  excite  both 
types  of  mode)  of  the  cavity  is  shown  in  Fig.  4b.  The  histogram  is  compared  with  two 
theories:  the  red  curve  corresponds  to  the  case  in  which  all  eigenfunctons  are  supposed  to 
be  of  the  chaotic  variety,  and  the  blue  curve  to  the  case  in  which  a  fraction  (as  determined 
by  theory)  are  of  the  integrable  type.  As  can  be  seen  the  theory  curve  including 
contributions  from  both  types  of  eigenfunction  matches  the  computed  distribution  of 
impedance  values.  We  are  currently  generalizing  the  formulation  of  this  effect  to  3D 
electromagnetic  cavities. 


Figure  4:  Modes  of  mushroom  shaped  cavity,  a)  two  adjacent  modes  showing  reflecting  chaotic  and 
regular  behavior  of  classical  ray  trajectories,  b)  Histograms  of  impedance  for  excitation  by  a  port  located  in 
the  mushroom  cap. 

Statistics  of  tunneling 

In  collaboration  with  L.  Pecora  of  NRL  we  have  developed  a  theory  for  the 
coupling  of  two  wave  enclosures  by  a  barrier  through  which  waves  must  tunnel.  [18]. 
The  main  conclusions  of  this  work  is  that  the  mean  tunneling  rate,  which  characterizes 
the  coupling  between  the  two  enclosures  is  independent  of  enclosure  shape  and  only 


depends  on  properties  of  the  tunneling  barrier.  However,  the  fluctuations  in  the  tunneling 
rate  from  mode  to  mode  depend  strongly  on  the  enclosure  shape.  Enclosures  with  chaotic 
ray  trajectories  have  much  smaller  fluctuations  than  enclosures  with  nonchaotic 
trajectories.  This  information  can  aid  in  design  or  placement  of  apertures  connecting 
enclosures  to  minimize  large  coupling.  This  general  theory  will  be  tested  experimentally 
in  the  proposed  work  as  discussed  below. 

Improving  TR  wave  focusing  using  iterative  techniques 

Spatiotemporal  focusing  of  waves  has  applications  in  fields  such  as  imaging  and 
communication,  but  is  also  extremely  useful  in  inducing  HPM  effects.  Time  reversal 
(TR)  mirrors  have  been  used  to  focus  waves  in  both  space  and  time  [19],  as  discussed 
above.  In  practice,  TR  mirrors  have  several  limitations  that  result  in  loss  of  information 
about  the  waves  broadcast  by  the  source;  these  include  i)  limited  coverage  by  the 
transceivers,  and  ii)  dissipation  during  the  wave  propagation  (which  breaks  TR 
invariance)  [20,  21].  The  limitation  due  to  dissipation  leads  to  increasing  loss  of 
information  as  the  recording  time  increases. 

The  loss  of  information  during  the  reconstruction  results  in  temporal  and  spatial  sidelobes 
of  the  reconstructed  pulse.  In  previous  work,  we  used  the  compensating  technique  of 
exponential-in-time  amplification  of  the  rebroadcasted  time-reversed  signal  to  partially 
undo  the  adverse  effects  of  dissipation,  and  to  enhance  the  range  of  sensors  which  utilize 
TR  mirrors  [20,  22].  However,  this  technique  does  not  improve  the  temporal  focusing  of 
the  reconstructed  pulse.  On  the  other  hand,  an  iterative  TR  technique  has  been  shown  to 
be  effective  in  eliminating  the  spatiotemporal  sidelobes  of  the  reconstructed  pulse  [23]. 
We  introduce  into  the  iteration  method  an  amplification  parameter  to  compensate  for 
dissipation.  By  tuning  this  parameter,  we  can  substantially  improve  the  accuracy  and 
convergence  of  the  iterative  focusing  technique.  This  is  demonstrated  both 
experimentally  and  numerically  [11]. 

Sensing  small  perturbations  using  TR  focusing 

We  demonstrated  a  remote  sensor  scheme  by  applying  the  quantum  mechanical 
concept  of  fidelity  loss  to  classical  acoustic  and  electromagnetic  waves  [8,  9].  The  sensor 
makes  explicit  use  of  time-reversal  invariance  and  spatial  reciprocity  in  a  wave  chaotic 
system  to  sensitively  and  remotely  measure  the  presence  of  small  perturbations.  The  loss 
of  fidelity  is  measured  by  employing  a  single-channel  time-reversal  mirror  to  rebroadcast 
a  probe  signal  into  the  perturbed  system.  We  also  introduced  the  use  of  exponential 
amplification  of  the  probe  signal  to  partially  overcome  the  effects  of  propagation  losses 
and  to  vary  the  sensitivity. 

We  then  showed  that  the  sensor  can  be  used  to  quantify  a  class  of  perturbations 
using  classical  waves.  In  particular,  we  investigated  perturbations,  which  change  the 
volume  of  a  strongly  scattering  wave  chaotic  cavity,  but  preserve  its  shape.  The  crux  of 
this  approach  is  that  the  effect  of  a  volume  changing  perturbation  is  a  change  in  the 
frequency  spacing  of  the  modes  of  the  cavity.  Thus  the  effect  of  the  perturbation  on  a 
response  signal  from  a  perturbed  cavity  can  be  undone  by  empirically  scaling  the  phase- 
information  contained  along  the  time-axis  of  the  response  signal  from  the  perturbed 
cavity.  The  scaling  value  that  undoes  the  perturbation,  and  hence  quantifies  the 
perturbation,  is  the  ratio  of  the  mode  spacing  before  and  after  the  perturbation.  The 
capability  to  quantify  such  volume  changing  perturbations  has  been  demonstrated  both 
numerically  and  experimentally.  A  quasi- ID  star  graph  model  [11]  was  initially  used  to 


demonstrate  the  results.  The  quantification  of  electrical-volume  changing  perturbations  to 
a  one  cubic  meter  aluminum  box  was  demonstrated  experimentally;  the  experimental 
results  are  also  supported  by  a  finite  difference  time  domain  simulation  of  the  full-wave- 
solution  to  the  time-domain  Maxwell's  Equations  inside  the  box.  Finally,  the  approach  to 
quantify  volume  changing  perturbations  was  shown  to  apply  to  a  generic  wave  chaotic 
system  by  using  a  time  domain  version  of  our  Random  Coupling  Model  [24] . 

Directing  wave  energy  using  TR  techniques 

Nonlinear  Time-Reversal  Mirror  -  An  ideal  time-reversal  mirror  in  an  open 
environment  would  collect  the  forward-propagating  wave  at  every  point  on  a  closed 
surface  enclosing  the  transmitter,  requiring  a  very  large  number  of  receivers.  The 
receiving  array  can  be  simplified,  without  significant  loss  of  fidelity  of  the  reconstruction, 
if  there  is  a  closed,  ray-chaotic  environment  where  a  propagating  wave  (with  wavelength 
Nonlinear  Time-Reversed  Pulse 

Collapses  at  the  Diode!  Sona  with  Nonlinearity 


Figure  5  Illustration  of  the  nonlinear  time  reversal  experiment  [22].  A  pulse  is  injected  at  one  port  of  an 
enclosure.  Signals  in  the  enclosure  interact  with  a  diode  and  are  then  collected  at  a  second  port.  When  that 
signal  is  time  reversed,  filtered,  and  reinjected,  a  pulse  will  form  cither  at  the  original  port  or  the  diode 
depending  on  the  filter. 

much  smaller  than  the  size  of  the  enclosure)  will  eventually  reach  every  point  in  the 
environment,  allowing  the  use  of  a  single  receiver  to  capture  the  signal  to  be  time- 
reversed.  Reconstruction  is  possible  even  when  only  a  small  fraction  of  the  transmitted 
energy  is  collected  by  the  receiver.  Now  consider  a  discrete  nonlinear  element  added  to 
the  otherwise  linear  environment.  When  a  waveform  is  incident  on  the  nonlinear  element, 
excitations  are  formed  at  frequencies  different  from  those  in  the  initial  pulse.  These  new 
excitations  appear  as  a  new  transmission  originating  from  the  nonlinear  element,  which  in 
principle  should  be  time-reversible  in  their  propagation,  and  behave  similarly  to  the 
initial  pulse.  In  particular,  the  nonlinear  excitations  should,  upon  time-reversal, 
reconstruct  as  a  signal  that  is  focused  on  the  nonlinear  element.  We  have  demonstrated 
this  method  of  nonlinear  time-reversal  with  electromagnetic  waves  [12]  as  illustrated  in 
Fig.  5.  Furthermore,  the  generation  and  time-reversal  of  nonlinear  excitations  will  not 
depend  upon  the  location  of  the  object,  which  may  be  unknown.  This  allows  creation  of 
an  exclusive  communication  channel  with  the  nonlinear  object,  and  the  ability  to  direct 
intense  energy  onto  it  without  interfering  with  nearby  objects. 

Experimental  studies  of  circuit  effects 

Under  previous  funding  a  comprehensive  program  of  effects  tests  on  electronic 
systems  lead  to  characterization  of  fundamental  effects  mechanisms  from  which 
generalized  effects  models  were  developed.  A  key  aspect  is  illustrated  in  Fig.  6,  which 


shows  a  typical  CMOS  integrated  circuit  along  with  a  comparison  of  results  from 
measurements  and  simulations.  HPM  coupled  onto  the  chip  is  likely  to  drive 
semiconductor  junctions  into  nonlinear,  non-quasistatic  (NQS)  regimes  of  operation.  We 
demonstrated  that  a  modification  to  the  Berkeley  (BSIM)  device  macro-models 
accurately  predicts  HPM  effects  in  electronic  circuits  using  commercially  available 
circuit  simulators  such  as  CADENCE  and  Agilent  ADS.  Our  approach  involves 
fabricating  basic  integrated  circuits  (IC)  and  measuring  their  response  when  driven  by 
continuous,  pulsed  and  wideband  HPM  waveforms.  The  results  allow  us  to  characterize 
the  various  linear  and  nonlinear  effects  transfer  functions  from  which  scalable  (device- 
based)  macro-models  are  constructed. 

We  have  shown  experimentally  that  HPM  fields  couple  predominantly  to  bus 
wires  and  board-level  interconnects  and  therefore  may  be  considered  as  radiation 
impedances  (e.g.  as  input  to  RCM).  Since  the  IC  dimensions  (mm's)  are  much  smaller 
than  typical  HPM  wavelengths  (10’s  cm),  the  IC  is  assumed  to  be  a  lumped  element 
excited  by  the  terminal  voltages  of  the  interconnects.  Everything  up  to  this  point  in  the 
model  is  linear.  Once  the  HPM-induced  voltage  levels  and  spectra  are  known  the 
nonlinear  device  models  discussed  above  are  employed  to  calculate  effects. 


(a)  (b) 

Figure  6:  Shown  is  (a)  the  VLSI  layout  of  a  CMOS  integrated  circuit  and  (b)  a  comparison  of  the  results 
of  measurements  and  simulation  that  include  a  modified  MOS  channel  model  that  correctly  predicts  NQS 
behavior  due  to  HPM  excitation. 


The  capabilities  described  above  will  now  enable  us  to  design,  evaluate  and  test 
new  concepts  for  HPM-hardened.  We  have  already  evaluated  the  HPM  susceptibility  of 
common  and  differential  mode  signaling  buses  and  found  that  differential  busses  are  1 7 
dB  more  robust  than  SE.  Similar  improvements  aimed  at  hardening  devices,  IC  boards, 
system  interconnects  and  electronic  enclosures  are  feasible.  Altogether  it  is  likely  that 
vastly  superior  HPM  immunity  in  mission-critical  electronic  systems  may  be  achieved. 
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This  thesis  treats  two  general  problem  areas  in  the  field  of  wave  chaos. 

The  first  problem  area  that  we  address  concerns  short  wavelength  tunneling 
from  a  classically  confined  region  in  which  the  classical  orbits  are  chaotic.  We  de¬ 
velop  a  quantitative  theory  for  the  statistics  of  energy  level  splittings  for  symmetric 
chaotic  wells  separated  by  a  tunneling  barrier.  Our  theory  is  based  on  the  ran¬ 
dom  plane  wave  hypothesis.  While  the  fluctuation  statistics  are  very  different  for 
chaotic  and  non- chaotic  well  dynamics,  we  show  that  the  mean  splittings  of  differ¬ 
ently  shaped  wells,  including  integrable  and  chaotic  wells,  are  the  same  if  their  well 
areas  and  barrier  parameters  are  the  same.  We  also  consider  the  case  of  tunneling 
from  a  single  well  into  a  region  with  outgoing  quantum  waves. 

Our  second  problem  area  concerns  the  statistical  properties  of  the  impedance 
matrix  (related  to  the  scattering  matrix)  describing  the  input/output  properties  of 
waves  in  cavities  in  which  ray  trajectories  that  are  regular  and  chaotic  coexist  (i.e., 
‘mixed’  systems).  The  impedance  can  be  written  as  a  summation  over  eigenmodes 


where  the  eigenmodes  can  typically  be  classified  as  either  regular  or  chaotic.  By 
appropriate  characterizations  of  regular  and  chaotic  contributions,  we  obtain  statis¬ 
tical  predictions  for  the  impedance.  We  then  test  these  predictions  by  comparison 
with  numerical  calculations  for  a  specific  cavity  shape,  obtaining  good  agreement. 
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Chapter  1 


Introduction 
1.1  Introduction 

According  to  the  correspondence  principle,  the  predictions  of  quantum  and 
classical  mechanics  should  coincide  in  the  limit  of  short  quantum  wavelength. 

It  is  particularly  interesting  to  investigate  possible  manifestations  of  the  corre¬ 
spondence  principal  in  situations  where  the  quantum  and  classical  pictures  display 
seemingly  different  fundamental  properties.  For  example,  classical  mechanics,  being 
nonlinear,  may  commonly  yield  chaos,  while  quantum  mechanics,  e.g.,  as  described 
by  the  Schrodinger  wave  equation,  is  linear  and  thus  cannot  yield  chaotic  dynamics 
in  the  usual  classical  sense  (exponential  sensitivity  of  bounded  solutions  to  small 
perturbation).  Thus,  the  short  wavelength  quantum  manifestations  of  chaos  in  a 
corresponding  classical  system  has  attracted  much  attention  [1,  2,  3,  4,  5],  and  the 
study  of  this  issue  has  been  given  the  appellation  ‘quantum  chaos’. 

An  early  consideration  that  later  turned  out  to  be  important  for  wave  chaos 
was  provided  by  Wigner  who  considered  energy  levels  of  large  nuclei  [6,  7,  8,  9,  1,  2], 
Since  the  energy  level  density  at  high  energy  is  rather  dense  and  the  solution  of 
the  wave  equation  for  the  levels  was  inaccessible,  Wigner  looked  for  a  statistical 
description  of  these  levels.  In  recent  years,  the  statistical  approach  to  properties  of 
the  solutions  of  wave  equations  in  systems  where  the  ray  trajectories  (e.g.,  classical 
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orbits  of  quantum  systems)  are  chaotic  has  been  a  very  active  area,  and  much  work 
has  been  done;  examples  include  optical  [10],  acoustic  [11],  microwave  [12,  13,  14,  15] 
and  electronic  cavities  [16,  17].  Here,  we  focus  on  quasi-two-dimensional  microwave 
cavities  and  quantum  dots,  both  of  which  are  assumed  to  be  thin  in  the  vertical 
direction  and  have  ray  trajectories  which  may  be  chaotic  in  the  other  two  dimensions 
(‘billiards’).  In  Chapter  2,  we  consider  what  happens  when  there  is  a  tunneling 
barrier  in  the  billiard  region,  while  in  Chapter  3  we  consider  coupling  to  an  external 
environment  through  suitable  openings  (called  ‘leads’  or  ‘ports’). 

With  respect  to  our  work  in  Chapter  3,  we  note  that  statistical  properties 
in  chaotic  cavities  with  external  connections  have  been  well  studied  using  various 
approaches,  e.g.,  the  ‘Poisson  Kernel’  [18,  19,  20,  21,  14]  or  the  ‘Random  Coupling 
Model’  (RCM)  [12,  13,  14,  15].  However,  in  general,  such  systems  may  have  not  only 
either  all  chaotic  or  all  regular  orbits,  but  also  typically  have  a  mixture  of  coexisting 
chaotic  and  regular  orbits.  We  called  such  systems  ‘mixed’.  Mixed  systems,  in 
spite  of  their  wide  occurance,  have  been  little  studied  in  the  previous  literature  on 
wave  chaos.  Extending  the  application  of  the  RCM  to  mixed  systems  is  the  goal  of 
Chapter  3. 

With  respect  to  our  work  in  Chapter  2,  we  note  that,  in  addition  to  chaos, 
the  quantum  phenomenon  of  tunneling  through  classically  forbidden  regions  of  phase 
space  presents  another  striking  difference  between  quantum  and  classical  mechanics. 
Much  past  work  examining  the  issue  of  tunneling  in  classically  chaotic  systems  has 
been  done  (e.g.,  Refs.  [22,  23,  24,  25,  26,  27,  28,  29,  30,  31,  32,  33,  34,  35,  36]  and 
references  therein).  For  example,  one  question  that  has  been  extensively  studied  is 
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what  happens  when  a  quantum  state  is  initially  localized  in  an  integrable  region  of 
classical  phase  space  and  tunnels  through  to  a  chaotic  region  [22],  which  is  often 
called  dynamical  tunneling.  In  contrast,  in  Chapter  2,  as  in  Refs.  [33,  34,  37],  we 
consider  the  problem  of  quantum  tunneling  from  a  chaotic  region  through  a  classical 
barrier  in  the  absence  of  an  integrable  region.  While  Refs.  [33,  34,  37]  treat  this 
problem  in  the  case  of  smoothly  varying  potentials  with  spatially  narrow  tunneling 
paths  (e.g.  ‘instantons’),  our  concern  will  be  the  two-dimensional  case  where  there 
are  piecewise- uniform  potentials,  long  barriers,  and  confining  hard  walls  [38]  1. 

An  important  point  is  that  smooth  Hamiltonians  (e.g.,  H  =  (p2/2m)  + 
R(q)  with  R(q)  smooth  function  of  q)  with  completely  chaotic  dynamics  treated 
in  [33,  34,  37]  typically  do  not  occur  [although  they  may  be  thought  to  approximate 
systems  where  regular  regions  occupy  a  small  fraction  of  the  allowed  phase  space  vol¬ 
ume],  On  the  other  hand,  completely  chaotic  phase  space  dynamics  does  occur  for 
billiard  systems  (zero  potential  regions  bounded  by  hard  walls,  Sec.  1.2).  Because  of 
the  flexibility  of  billiard  systems  in  allowing  various  types  of  dynamics  (chaotic,  non- 
chaotic,  or  mixed)  this  thesis  will  concentrate  on  such  systems.  Other  motivations 
for  considering  billiard-type  systems  include:  (i)  they  are  potentially  realizable  in 
quantum  dot  contexts  and  in  descriptions  of  classical  optical  electromagnetic  holds 

in  piecewise-constant  dielectrics;  (ii)  by  adjusting  the  shape  of  the  billiards,  it  is 
1In  contrast  with  the  billiard- type  classical  chaos  that  we  treat,  smooth  Hamiltonians  typically 
yield  mixed  phase  spaces  with  coexisting  chaotic  and  regular  regions.  In  fact,  we  are  aware  of  only 
one  instance  where  a  smooth  Hamiltonian  system  has  been  claimed  to  be  fully  chaotic,  namely, 
the  anisotropic  Kepler  problem  with  zero  angular  momentum. 


3 


particularly  easy  to  go  from  integrable  to  mixed  to  chaotic  phase  space;  (iii)  com¬ 
parisons  between  different  shape  billiards  (e.g.,  between  an  integrable  and  a  chaotic 
case)  can  be  made  quantitatively  precise  by  keeping  certain  gross  parameters  equal 
(see  Sec.  2.1). 

For  example,  with  respect  to  point  (iii)  above,  a  major  result  of  Chapter  2 
is  that  integrable  and  chaotic  cases  with  the  same  mean  wave  tunneling  properties 
differ  very  greatly  in  their  fluctuation  characteristics,  with  the  chaotic  case  having 
much  smaller  fluctuations  about  the  common  mean  than  the  integrable  case.  We 
believe  that,  in  the  billiard  case  discussed  here,  point  (iii)  makes  this  effect  a  partic¬ 
ularly  dramatic  instance  of  a  quantum  manifestation  of  classical  chaos.  We  remark 
that  this  relative  suppression  of  tunneling  fluctuations  in  the  chaotic  case  occurs 
because,  due  to  the  classical  ergodicity  of  chaotic  systems,  the  quantum  states  are 
relatively  similar  in  that  they  typically  effectively  spread  over  the  entire  classically 
allowed  phase  space.  In  contrast,  in  integrable  systems,  due  to  the  existence  of 
extra  constants  of  the  motion,  different  energy  states  are  typically  constrained  to 
have  much  more  variation  of  their  distribution  in  phase  space  and  may  avoid  the 
phase  space  region  where  tunneling  is  strongly  excited.  If  so,  the  tunneling  can 
be  exponentially  small  and  very  dependent  on  the  particular  state.  This  point,  al¬ 
ready  inherent  in  the  discussions  in  Refs.  [33,  34,  37],  applies  to  both  the  case  of 
billiard  type  Hamiltonians  and  the  case  of  smooth  Hamiltonians  (e.g.,  see  Fig.  2  of 
Ref.  [33]). 
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(a)  (b) 


Figure  1.1:  First  100  bounces  of  two  trajectories  (blue  and  red)  started  at  the 
same  point  with  slightly  different  initial  directions  in  integrable  (a)square  billiard 
(b) quarter  circle  billiard. 

1.2  Billiard 


Classically,  a  billiard  is  a  dynamical  system  in  which  a  particle  moves  in  a 
straight  line,  with  constant  energy,  in  a  confined  region,  B,  and  is  reflected  specularly 
at  the  boundary,  dfl,  without  change  of  speed.  Billiard  systems  have  a  simple 
Hamiltonian, 

tf(p,q)  =  -||j  +  r(q),  (i.i) 


where 


10  for  q  G  B, 

(1.2) 

oo  for  q  d  B, 

but  the  particle  trajectory  inside  different  billiards  can  have  three  different  behav¬ 
iors  (a)  regular,  (b)  chaotic,  and  (c)  a  mixture  of  regular  and  chaotic.  The  shape 
of  the  boundary  determines  which  types  of  trajectories  exist  in  the  billiard  (see 
Figs.  1.1,  1.2,  and  1.3). 
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(a)  (b) 


Figure  1.2:  First  100  bounces  of  two  trajectories  (blue  and  red)  started  at  the  same 
point  with  slightly  different  initial  directions  in  chaotic  billiards  (a)  Sinai  billiard 
(b)  Stadium  billiard. 


Figure  1.3:  First  100  bounces  of  two  trajectories  (blue  and  red)  started  at  the  same 
point  with  slightly  different  initial  directions  in  the  mushroom  billiard  (a) chaotic 
(b)integrable. 
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Considering  the  quantum  billiard,  the  Hamiltonian  of  the  time-independent 
Schrodinger  equation,  Hcj)  =  E(f> ,  is  replaced  by  the  classical  Hamiltonian,  Eqs.  (1.1, 
1.2),  he., 

-;r-V20n(q)  =  En<t>n{ q)  for  q  e  H,  (1.3) 

2  m 

where  (j)n  and  En,  n  —  1,  2, . . . ,  oo  are  real  eigenfunctions  and  eigenvalues,  and  the 
infinite  potential  outside  the  region  translates  to  the  Dirichlet  boundary  conditions: 

<t>n( q)  =  0  for  q  £  fl  (1.4) 

The  eigenfunctions  are  chosen  be  orthogonal: 

[  0m(q)</>n(q)rfq  =  (i-5) 

Jn 

With  k2  =  2 mEnj h2,  this  free-held  Schrodinger  equation  is  the  same  as  the  Helmholtz 
equation, 

(V2  +  k2n )  M r)  =  0.  (1.6) 

The  relation  between  the  Schrodinger  equation  and  the  Helmholtz  equation 
implies  that  quantum  billiards  can  be  modeled  by  microwave  cavity  of  a  given  shape, 
thus  opening  a  door  to  experimental  verification.  In  particular,  we  note  that  the 
vertical  electric  held  of  modes  in  a  microwave  cavity  that  is  vertically  thinner  than  a 
half- wavelength  satisfy  (1.6)  in  two  dimensions  with  Dirichlet  boundary  conditions. 
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1.3  Spectrum 


1.3.1  Weyl’s  Formula  and  Normalized  Spacing 


For  a  d-dimensional  wave  equation,  (V2  +  A2) <f>  =  0,  in  a  region  £1  of  volume 
V,  the  number  N  of  eigenmodes  with  A2  below  k 2  can  be  described  by  the  Weyl’s 
formula  [39], 

N(k2)  =  (2  n)~dVkd  +  0{kd~l).  (1.7) 


Note  that  the  boundary  condition  of  the  Helmholtz  equation  will  have  a  small  effect 
of  0{kd~l)  in  the  mode  counting  formula  (1.7).  Specifically,  for  2-dimensional  wave 
equation  in  a  region  fl  of  area  A,  the  number  of  eigenmodes  below  than  k2  is 


N(k2)  =  A k 2  +  0(k), 


(1.8) 


or  the  mode  density  is 


P(k 2)  = 


A 
47 r 


(1.9) 


Thus,  the  average  spacing  between  eigenvalues  [A (k2)  =  k2+1  —  k2  for  k2  =  k2]  is 

1  47T 


A (k2)  = 


(1.10) 


p(k2)  A ' 

We  dehne  the  normalized  eigenvalue  spacing  using  the  average  spacing  between 


eigenvalues,  as  follows 


k2  —  k2 

_  '"'n+l  ft'n 

Sn  =  £ 


(1-11) 


1.3.2  Random  Matrix  Theory 

Random  matrices  were  introduced  by  Eugene  Wigner  to  model  the  spectra  of 
heavy  nuclei.  Since  the  energy  levels  at  high  energy  are  rather  dense,  and  the  wave 


equation  was  difficult  to  solve,  Wigner  sought  a  statistical  description  for  properties 
of  the  spectrum.  Wigner  hypothesized  that  the  eigenvalue  spectrum  of  complicated 
nuclei  have  similar  statistical  properties  to  those  of  the  spectra  of  ensembles  of 
random  matrices  that  depend  only  on  the  symmetry  of  the  Hamiltonian. 

This  led  Wigner  to  consider  three  types  of  random  matrix  ensembles,  called 
Gaussian  Orthogonal  Ensembles  (GOE),  Gaussian  Unitary  Ensembles  (GUE)  and 
Gaussian  Symplectic  Ensembles  (GSE).  In  this  thesis,  we  will  be  concerned  with 
GOE  ensembles  which  are  defined  to  have  the  following  two  properties. 

•  The  ensemble  is  invariant  under  every  orthogonal  transformation 

H  ->  OtHO  (1.12) 

where  O  is  any  real  orthogonal  matrix,  H7  =  H  and  H  is  real. 

•  The  various  elements  Hkj,  k  <  j,  are  statistically  independent  and 

{1V(0, 1)  for  i  —  j 

(1.13) 

N(  0,1/2)  for//./, 

where  N(n,a2)  is  a  Gaussian  random  variable  with  mean  /i  and  variance  a2. 

For  the  GOE  ensemble,  Wigner  found  that  the  distribution  of  normalized 
spacing  (see  Eq.  (1.11))  obeyed 

Pgoe(s)  =  |sexp  (-;|s2)  (1-14) 

(1.15) 

In  a  foundational  paper  for  quantum  chaos,  the  Bohigas-Giannoni-Schmit 
(BGS)  [40]  conjectured  that  the  short  wavelength  spectral  statistics  of  quantum 
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systems  whose  classical  counterparts  exhibit  chaotic  behavior  are  described  by  ran¬ 
dom  matrix  theory,  and  they  numerically  tested  that  (1.14)  was  satisfied  for  the 
billiard  of  Fig.  1.2. 

1.3.3  Level  Spacings  for  Regular  Systems 

In  a  classically  integrable  systems,  Berry  and  Tabor  [41]  showed  that  corre¬ 
sponding  solution  of  the  Schrodinger  equation  has  energy  levels  that  are  uncorrelated 
and  that  the  normalized  level  spacing  distribution  is 

Ppoisson^yS)  6  .  (1.16) 

1.3.4  Level  Spacing  for  Mixed  Systems 

Consider  a  system  whose  classical  phase  space  has  both  regular  and  chaotic 
regions,  and  denote  the  total  volume  ratio  of  the  regular  regions  by  pr  and  the 
volume  ratio  of  the  chaotic  region  by  pc.  Percival  conjectured  that,  in  the  semi- 
classical  limit,  the  energy  levels  consist  of  regular  and  chaotic  parts  having  certain 
distinct  properties  [42],  Berry  and  Robnik  extended  this  conjecture  and  assumed 
that  the  sequence  of  the  energy  levels  of  a  mixed  system  is  given  by  the  superposition 
of  statistically  independent  sub-sequences  corresponding  to  the  classical  phase-space 
regions  [43].  In  addition,  they  assumed  that  the  distributions  of  the  sub-sequences 
corresponding  to  regular  and  irregular  regions  are,  respectively,  the  Poisson  and 
Wigner  distributions,  see  Eqs.  (1.14)  and  (1.16). 
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1.4  Random  Plane  Wave  Hypothesis 


The  motivation  for  the  random  plane  wave  hypothesis  is  the  observation  that 
ray  trajectories  in  chaotic  billiards  (see  Figs.  1.2)  are  distributed  uniformly  in  space 
and  isotropicly  in  direction.  Correspondingly,  Berry  [44]  proposed  that  at  any  point, 
not  too  close  to  the  boundary,  the  eigenfunction  of  (V2  +  k^)(f)n  =  0  has  statistical 
properties  approximately  similar  to  those  of  a  random  superposition  of  many  plane 
waves  with  same  wavenumber,  kn, 

N 

0n(x)  =  otj  exp  ( iknej  ■  x  +  iOj )  +  (complex  conjugate),  IV  3>  1,  (1.17) 

3= 1 

where  the  amplitude  cp’s  are  identical  independently  distributed  random  variable 
with  some  probability  density  function,  the  direction  ej  are  independent  isotropicly 
distributed  random  unit  vectors,  the  phase  6j  are  independent  uniformly  distributed 
in  [0,  27t)  random  variables,  and  (1.17)  is  assumed  to  hold  when  2k /kn  is  small 
compared  to  a  typical  length  dimension  of  the  billiard.  Based  on  the  central  limit 
theorem,  it  is  thus  expected  that  the  eigenfunction  amplitude  0n(x)  at  any  given 
point  x  is  Gaussian  distributed  with  zero  mean  and  the  variance  is  1/V,  where  V 
is  the  volume  of  the  closed  region,  i.e., 


PM  = 


2i t<j2(x)  ^  [  2<v2(x)J  ' 

where  <r2(x)  =  1/V  (  u2(x)  =  1/V  follows  from  the  normalization  condition  (1.5)). 

From  the  random  plane  wave  hypothesis,  Berry  [44]  show  that  the  two-point 
correlation  function  in  a  d-dimensional  billiard  is 


(1.18) 


C(ri,r2',k)  =  , 


(1.19) 
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where  L  —  |rx  —  r2|,  V  is  the  d-dimensional  volume  of  the  billiard,  V  is  the  gamma 
function,  and  Jn  is  the  n-tli  order  Bessel  function  of  the  hrst  kind. 

We  compare  our  numerical  calculations  of  eigenfunctions  in  different  shape 
billiards  and  Berry’s  theory  in  Appendices  C  and  D. 

1.5  Outline  of  Dissertation 

This  dissertation  is  organized  as  follows:  In  Chapter  22,  we  study  the  statistics 
of  the  energy  level  splittings  between  symmetric  and  antisymmetric  pairs  of  mirror 
symmetric  wells  coupled  by  a  rectangular  tunneling  barrier.  We  use  the  random 
plane  wave  hypothesis  to  develop  a  theory  for  the  chaotic  cases.  We  also  show  that 
the  mean  splittings  of  differently  shaped  wells,  including  both  integrable  and  chaotic 
wells,  are  the  same  if  their  well  areas  and  barrier  parameters  are  the  same,  but  that 
the  statistics  of  fluctuation  are  very  different  for  chaotic  and  integrable  wells. 

In  Chapter  33,  we  study  the  statistics  of  the  input/output  properties  of  waves 
in  mixed  cavities  in  which  ray  trajectories  that  are  regular  and  chaotic  coexist.  In 
particular,  we  focus  on  the  statistical  properties  of  the  impedance  matrix  (related 
to  the  scattering  matrix)  which  can  be  written  as  a  sum  over  eigenmodes  where  the 
eigenmodes  can  typically  be  classified  as  either  regular  or  chaotic.  We  obtain  statis¬ 
tical  predictions  for  the  impedance  by  separating  the  regular  and  chaotic  contribu- 
2Chapter  2  is  a  republication  of  work  published  in  Physical  Review  E,  as  approved  by  the  thesis 
committee  [45]. 

3Chapter  3  is  a  republication  of  work  published  in  Physical  Review  E,  as  approved  by  the  thesis 
committee  [46]. 
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tions.  Finally,  we  test  these  predictions  by  comparison  with  numerical  calculations 
for  a  specific  mushroom  cavity  shape  and  obtain  good  agreement. 
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Chapter  2 


Theory  of  Chaos  Regularization  of  Tunneling  in  Chaotic  Quantum 
Dots 

2.1  Introduction 

The  work  presented  in  this  chapter  is  a  continuation  of  previous  work  in 
Ref.  [47]  in  which  we  reported  numerical  results  and  abbreviated  heuristic  argu¬ 
ments  justifying  our  numerical  observations.  Our  aim  now  is  to  provide  a  fuller 
theoretical  analysis  of  the  results  in  Ref.  [47]. 

Reference  [47]  considers  symmetric  double  well  situations  of  the  type  shown  in 
Figs.  2.1  (a)  and  (b).  There  is  a  barrier  region  of  uniform  potential  Vb,  width  2A, 
and  length  L.  This  barrier  region  separates  two  mirror-symmetric  wells  in  which  the 
potential  is  zero  and  whose  (non-barrier)  boundaries  are  hard  walls.  If  the  energy  E 
is  less  than  Vb,  then  a  point  particle  is  classically  confined  to  one  of  the  wells,  and 
its  orbit  follows  straight  lines  between  specular  reflections  from  the  well  boundaries 
(a  billiard  system).  The  character  of  the  orbit  depends  on  the  shape  of  the  well. 
For  the  rectangular  well  of  Fig.  2.1  (a)  the  orbits  of  a  point  particle  are  integrable, 
with  separately  constant  horizontal  and  vertical  kinetic  energies.  For  the  shape  of 
the  well  in  Fig.  2.1  (b),  the  convex  walls  insure  that  all  typical  orbits  are  chaotic 
and  ergodically  fill  the  full  available  phase  space  [47].  In  particular,  if  a  typical 
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INTEGRABLE  CHAOTIC 


Figure  2.1:  Symmetrical  double  wells  of  area  A  separated  by  a  tunneling  potential 
barrier  of  width  2  A.  length  L  and  height  Vb ■  (a)  shows  the  case  of  rectangular  wells, 
while  (b)  shows  a  case  in  which  all  typical  orbits  are  chaotic. 

particle  orbit  in  the  Fig.  2.1  (b)  billiard  is  sampled  at  some  random  time  t,  then  the 
location  of  the  particle  has  a  uniform  probability  density  per  unit  area  in  the  well, 
and  the  probability  density  of  the  direction  of  particle  motion  is  uniformly  isotropic 
in  [0,  27t).  (Reference  [47]  also  treats  other  completely  regular  or  chaotic  well  shapes 

l 

Considering  symmetric  wells,  as  in  Fig.  2.1,  the  quantum  eigenstates  have  ci¬ 
ther  even  or  odd  parity  with  respect  to  the  center  line,  and  for  E  sufficiently  below 
the  barrier  height  Vb,  we  may  consider  the  states  to  come  in  symmetric/antisym¬ 
metric  pairs  with  nearly  equal  energies.  We  denote  the  o  th  such  pair  (E^,E^). 

The  symmetric  state  energy  is  always  less  than  the  antisymmetric  state  energy, 
1These  include  the  stadium  billiard,  which,  although  chaotic,  due  to  its  continuous  family  of 
neutrally  stable,  ’bouncing-ball'  periodic  orbits,  exhibits  scar-type  modes  with  tunneling  rates  that 
can  substantially  deviate  from  the  mean.  Note  that,  due  to  its  convex  walls,  such  orbits  are  absent 
in  the  strongly  chaotic  case  of  Fig.  2.1(b) 
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Figure  2.2:  Energy  level  splittings  versus  energy  plotted  as  black  dots,  along  with  the 
sliding  average  (red),  (A E)E)t.  The  parameters  used  in  there  plots  are  Vb  =  1000, 
Ab  =  0.05,  L  =  2.423,  A  =  4.8. 

E s  <  E The  energy  level  splitting  is  denoted 

A E„  =  E, A  -  El  (2.1) 

Figure  2.2  shows  as  black  dots  values  of  A Ea  versus  Ea  =  \{E^+E^)  from  numerical 
solutions  of  the  normalized  Schrodinger  equation, 

{\72  +  E -V(x,y)}i>(x,y)  =  0,  (2.2) 

with  ^  =  0  on  the  hard  walls,  V  —  Vb  in  the  barrier  region  (0  <  x  <  2A)  and  V  =  0 
in  the  wells.  The  parameters  Vb,  A  and  L  and  the  well  area  A,  are  all  taken  to  be 
the  same  for  the  two  cases,  Fig.  2.1  (a)  and  Fig.  2.1  (b).  Also  EA  1,  i.e.,  the 


16 


well  dimensions  are  large  compared  to  the  quantum  wavelength,  corresponding  to 
the  semiclassical  regime.  Shown  in  Fig.  2.2  (a)  and  (b),  in  red,  is  a  sliding  average 
(A E)e,s  of  A Ea  using  a  window,  ( E  —  e)  to  (E  +  e),  that  encompasses  2  to  15 
splitting  values.  Figure  2.3  shows  the  two  sliding  averages  plotted  together  on  the 
same  graph,  in  blue  for  the  integrable  case  (Fig.  2.1  (a))  and  in  black  for  the  chaotic 
case  (Fig.  2.1  (b)),  along  with  our  theoretical  result  (red)  to  be  derived  in  Sec.  2.5. 
The  two  main  conclusions  from  the  numerical  results  of  Ref.  [47]  are  the  following: 

(1)  Fluctuations  of  the  quantum  splittings  are  very  much  larger  (note  the 
logarithmic  vertical  scale)  in  the  integrable  well  case  as  compared  to  the  chaotic 
well  case. 

(2)  For  the  same  gross  parameters  (A,Vb,  A,  L)\  the  sliding  average  (A E)e,e 
versus  E  is  independent  of  the  well  shape. 

In  Ref.  [47]  it  was  found  that  (1)  and  (2)  hold  for  all  pairs  of  similarly  related 
chaotic  and  regular  well  shapes  studied  1 . 

The  rest  of  this  chapter  is  organized  as  follows.  Noting  that  the  numerical 
results  in  Fig.  2.2  all  satisfy 

A Ea  <C  -  (Ea+ 1  —  Ea_ i) ,  (2-3) 

in  Sec.  2.2  we  use  perturbation  theory  to  develop  a  formal  expression  for  A Ea. 
This  expression  is  essentially  Herring’s  formula  [48].  Herring’s  formula  was  also 
used  in  Wilkinson’s  treatment  of  tunneling  between  regular  regions  adapted  to  our 
problem  [33].  Section  2.3  then  applies  Berry’s  random  plane  wave  hypothesis  [44] 
to  obtain  the  statistics  of  the  splittings  {AEa}  in  the  case  of  chaotic  classical  dy- 
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E 

Figure  2.3:  Sliding  average  versus  E. 


namics.  In  Sec.  2.4,  as  an  example,  we  numerically  test  the  result  of  Sec.  2.3  by 
comparing  its  predictions  with  the  data  shown  in  Fig.  2.2  (b).  Section  2.5  applies 
a  Green’s  function  technique  based  on  the  method  of  Balian  and  Bloch  for  deriving 
the  semiclassical  perimeter  correlation  to  the  density  of  state  for  billiard  systems 
[49,  50,  51,  52]  to  obtain  the  sliding  average  splitting  (A E)Eje  and  show  that  it  is  in¬ 
dependent  of  well  shape.  Section  2.6  concludes  with  further  discussion.  As  discussed 
in  Sec.  2.6,  we  report  in  Sec.  2.6.2  the  extension  of  our  method  to  the  treatment  of 
tunneling  out  of  a  single  chaotic  well  into  an  open  region  with  outgoing  quantum 
waves. 
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2.2  Perturbation  Theory  for  the  Statistics  of  Energy  Level  Splitting 


for  Chaotic  Wells 


2.2.1  Setup 

We  consider  the  symmetric  and  antisymmetric  wave  functions,  denoted  ipg 
and  tpA,  along  with  their  corresponding  energy  levels,  Eg  =  k2s  and  E2 4  =  kA  (where 
we  choose  units  in  which  K2 /(2m)  =  1).  Referring  to  Fig.2.1,  we  take  the  potential 
to  be  zero  in  the  left  and  right  wells  and  to  be  V  =  Vg  =  k2B  in  the  barrier  region, 
0  <  x  <  2A.  Focusing  on  the  left  well,  we  have  that 


and 


'd_ 

dx 


(V2  +  k%A)il> 
P sAxiV ) 


g  A  =  0  for  x  <  0, 

+  HsA^sA^  ,y)\  =  0, 


rr=0_ 


(2.4) 

(2.5) 


*PsAx,y)  =  °  (2.6) 

on  the  boundary  of  the  left  well  other  than  that  at  x  =  0.  Hg  and  HA  denote 
operators  on  functions  of  y  that  we  now  obtain. 


2.2.2  Boundary  Condition  at  x  =  0 

Within  the  barrier  region,  0  <  x  <  2A,  -0a  and  bs  satisfy  the  following 
conditions, 

-  =  0  and  ip  a.  =  0  at  x  —  A.  (2.7) 

ox 
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Thus,  in  the  barrier,  solutions  of  the  time-independent  Schrodinger  equation, 


[VL  -  (krB  -  k2S)A)}^s,A  =  o,  i>s,A  =  0  at  y  =  0,  L, 


can  be  written  as 


^s{x,y)  =  ^2  Sr 


cosh  [am(A  —  &)]  .  f  rmry\ 


m=  1 
oo 


cosh  (am A) 


sm 


V  L  )' 


,  ,  \  a  sinh  k(A  - ^)]  •  /'miry 

ipA{x,  y)  =  2_ A ™ - — —  sm 


m=l 


sinh  (am A) 


V  L 


Oim  =  \lk%  +  (r^-'\  -k2. 


V  L  ) 


Noting  that  both  i/)a,s  and  d^A,s/dx  are  continuous  at  x  =  0 
Hs  and  Ha  in  (2.5)  are  given  by 

rv~\ 

r(m)  f  0-  ( rrmy 


Hs[f{y )]  =  Hg  fm  sin  [  L  n 

m=  1 

Hgn)  =  am  tanh  (amA), 

OO 

HA[f(y)]  =  ^(rrvK^- 


sm 


m—  1 


H =  am  coth(amA). 


where  fm  denote  coefficients  of  the  Fourier  sine  series  of  f(y), 


f(y)  =  fm  sin  ( 


/  rrmy 


m=  1 


Note  also  that  the  H,s,a  operators  are  self  adjoint, 

[  g(y)HsAf(y))dy  =  f  f(y)Hs,A\g(y)]dy. 
Jo  Jo 


(2.8) 

(2.9) 
(2.10) 
(2.11) 
,  we  have  that 

(2.12) 

(2.13) 

(2.14) 

(2.15) 
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2.2.3  Perturbation  expansion 


As  the  thickness  of  the  barrier  A  becomes  large,  we  see  from  (2.12)  and  (2.13) 
that  Hg  and  Ha  become  the  same: 

f)W  fr(m) 

12  S  '  !I  A  f 

We  denote  this  limit  by  the  subscript  0  and  define  a  corresponding  wavefunction 
and  energy  level,  bo  and  k q,  that  satisfy  the  problem, 

(V2  +  /cg)b o  =  0  for  x  <  0,  (2.16) 

+  H0[MO~,y)]  =0,  (2.17) 

x=0~ 

and  bo  =  0  on  the  boundaries  of  the  left  well  other  than  that  at  x  —  0.  The  operator 
Hq  is  defined  as  in  Eqs.  (2.12)  and  (2.13)  with 

=  am.  (2.18) 

Since  Hg  and  Ha  become  equal  as  A  — >  oo,  the  symmetric  and  antisymmetric  energy 
eigenfunctions  (bs  and  b^  )  and  energy  levels  ( k2s  and  k a)  also  become  equal.  (In 
particular,  they  become  bo  and  k g.)  Thus,  for  sufficiently  large  A,  we  can  assume 
that  these  symmetric  and  antisymmetric  quantities  are  close  to  each  other  and  are 
close  to  the  solution  of  Eqs.  (2.16)-  (2.18).  More  formally,  if  we  introduce  a  small 
expansion  parameter  e,  we  have 

b s,a  -  bo  =  0(e) 
kl,A  ~kl  =  0(e) 

HSA-H0  =  O(e).  (2.19) 
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Multiplying  Eq.  (2.4)  by  'ip0(x,y)dxdy  and  Eq.  (2.16)  by  tps,A(x,y)dxdy,  integrating 
over  the  area  of  the  left  well,  and  subtracting  the  results,  we  obtain 


/  [0O V20s,a  -  0s, a V20o]  dxdy  =  {k2SA  -  k20)  /  0O ^s^dxdy,  (2.20) 

Jlw  ’  Jlw 

where  LW  denotes  the  area  of  the  left  well.  Applying  Green’s  identity  to  the  left 

side  of  this  equation,  we  have  essentially  Herring’s  formula  [48], 


{-00 Hs,a[^s,a\  ~  0s,A#o[0o]}*=och/  =  {k%A  ~ 


0o  i>s,Adxdy,  (2.21) 


J  o  Jlw 

where  we  have  used  the  condition  that  ips,Afl  —  0  on  the  boundaries  of  the  left  well 
other  than  that  at  x  —  0.  Furthermore,  from  (2.15)  the  left  side  of  (2.21)  can  be 
rewritten 


/  {;i})^Hs.a bks.A] }x=ody  =  {k2S  A  -  kl )  /  0o i/js,Adxdy,  (2.22) 

Jo  Jlw 

where  A Hs,A[f(y)]  =  HsAfiv )]  “  H0[f(y)].  Noting  from  (2.19)  that  A HSA  and 

{k2s A  —  kl)  are  both  O(e),  we  see  that,  to  lowest  order  in  e,  we  can  set  ips,A  =  0o 

in  (2.22),  thus  yielding  the  perturbation  theory  result, 


,2  1.2  _  Jo  {00 [0o]  }x=ody 

S’A  °“  Lw^dxdy  ■ 


(2.23) 


ft  follows  from  (2.12),  (2.13)  and  (2.18)  that  >  H^1'1  >  Hg,l\  Thus,  we  have 
that  kA>  kl  >  kg.  We  denote  the  energy  level  splitting  by  A E  =  A k2  =  k\  —  k2s. 
Taking  the  difference  between  the  symmetric  and  antisymmetric  versions  of  (2.23) 
and  employing  our  definitions  of  JJs,a, o,  we  obtain 


A  E  =  LamPln  ,  f 

“  sinh  (2amA) 7  JLW 


i^dxdy, 


(2.24) 
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where  Cm  are  the  Fourier  sine  coefficients  of  i/)o(0,y), 


V’o(O)  y)  =  X]  C'm  sin  ■  (2-25) 

m=l 


2.3  Evaluation  of  A E  for  Chaotic  Eigenfunctions 


In  the  region  x  <  0,  but  near  x  —  0,  the  upper  well  boundary  in  Fig.  2.1  is 
close  to  y  =  L.  Therefore,  in  this  region  we  can  take 

OO  _ 

V>oO,  y)  =  ^2  cm  sin  sin  (kx.mx  -  (2.26) 

m=l 

where  k^m  =  k$  —  ( rmr/L )2,  and  4)rn  is  determined  by  the  boundary  condition 
d'l'pQ /dx  =  —  Hofyol  at  x  —  0,  which  yields 

tan0m  =  ^2h.  (2.27) 

Comparing  (2.26)  and  (2.25),  we  see  that 


Crn  —  — Cm  sill  ( 


(2.28) 


and  from  (2.11)  and  (2.27) 

sin0m  =  —j—.  (2.29) 

ks 

We  will  view  Eq.  (2.26)  as  is  a  statistical  model  and  think  of  the  values  of  the 
cm  as  pseudo-random  variables  that,  for  any  given  two  realizations,  can  be  regarded 
as  representing  two  different  eigenfunctions  with  nearly  the  same  energy  k'fy  In  what 
follows  we  will  approximate  (2.26)  by  cutting  off  the  sum  at  m  =  M,  where  M  is 
the  maximum  value  of  m  for  which  k%  m  >  0, 


M 


" 

,  .  /m7TM 

max 

m 

Al 

O 

■se 

(2.30) 
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That  is,  we  only  include  propagating  modes. 

We  now  need  a  model  for  characterizing  the  pseudo-random  coefficients  cm 
in  (2.26).  To  do  this,  we  assume  M  3>  1,  follow  Berry  [44],  and  utilize  the  chaotic 
classical  dynamics  of  particles  in  the  potential  wells,  together  with  the  correspon¬ 
dence  principle.  Our  chaotic  classical  particle  trajectories  have  the  following  ergodic 
character:  For  typical  initial  conditions  and  any  small  localized  region  SR  in  the  well, 
a  very  long  orbit  will  pass  through  SR  many  times,  and,  if  one  examines  these  pas¬ 
sages,  one  will  find  that,  as  the  orbit  length  approaches  infinity,  (i)  the  fraction  of 
time  spent  by  the  orbit  in  the  region  SR  is  the  ratio  of  the  area  of  SR  to  the  total  area 
of  the  well,  and  (ii)  the  orientation  of  the  particle’s  velocity,  in  its  passes  through 
SR,  is  uniformly  distributed  in  angle.  Thus,  if  we  imagine  sampling  the  chaotic  orbit 
at  some  randomly  chosen  time,  its  location  will  have  a  uniform  probability  density 
distribution  in  space  and  its  velocity  (whose  magnitude  is  fixed  by  the  particle  en¬ 
ergy)  will  have  an  isotropic  probability  distribution  in  its  orientation.  Thinking  of 
-00  as  analogous  to  the  classical  probability  density  in  space  and  invoking  property 
(i),  the  correspondence  principle  suggests  that,  for  wavelengths  small  compared  to 
the  cavity  size,  the  coarse  grained  average  of  over  several  wavelengths  will  have 
a  value  that  is  the  same  near  x  =  0  as  in  any  other  region  of  the  well. 

We  now  ask  how  the  coefficients  cm  in  (2.26)  are  related  to  the  integral 
fLw^dxdy  appearing  in  the  denominator  of  Eq.  (2.24).  We  see  from  (2.26)  that 
the  integral  of  0q  over  a  region  of  area  Ab  abutting  the  barrier  and  extending  not 
too  far  away  from  it  is  {Ab/ 4)  where  the  factor  1/4  results  from  taking  the 

spatial  average  of  sin2  {rrmy/L)  sin2  (kx,mx  —  <f>m)  over  several  wavelengths.  If  the 
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wavelength  is  short,  k0L  S>  1,  based  on  point  (i)  above  and  the  correspondence 
principle,  one  might  suppose  that  this  result  for  fA  -0g dxdy  can  be  extended  to  x 
values  far  from  the  barrier.  In  particular,  we  expect  that  the  spatial  average  of  ?/;,j 
near  x  =  0  is  approximately  the  same  as  in  any  other  region  of  the  well.  Thus,  we 
obtain  the  estimate, 


[  =  (2-31) 

J  LW  ^ 

vy  m 

where  (as  previously  stated)  A  is  the  entire  area  of  the  left  or  right  well.  We 
expect  (2.31)  to  hold  as  long  as  the  barrier  dimension  L  is  much  greater  than  a 
wavelength  M 7r  ~  k^L  3>  1.  When  the  barrier  dimension  becomes  comparable  to 
or  smaller  than  a  wavelength,  point  to  point  variations  in  the  magnitude  of  the 
wavefunction  lead  to  departures  between  the  average  value  of  in  the  well  and 
the  corresponding  value  near  the  barrier.  However,  even  if  k$L  is  large,  Eq.  (2.31) 
is  only  approximate,  and  we  expect  it  to  hold  with  an  error  that  becomes  small  as 
k0L  — >  oo.  We  will  find  that,  when  computing  fluctuations  in  energy  splittings,  the 
small  fluctuating  error  in  (2.31)  can  be  important  (e.g.,  see  Sec.  2.6.1  and  Sec.  2.4). 
Use  of  (2.31)  as  a  strict  equality  assumes  that  the  coarse  grained  average  of  i/jq  is 
essentially  determined  throughout  the  total  area,  A,  by  the  M  amplitudes,  Cm  of 
the  propagating  modes  (k%  m  >  0)  near  the  barrier.  This  is  not  always  the  case.  For 
example,  if  scars  are  present,  there  may  be  deviation  between  the  average  of  ^  near 
the  barrier  and  throughout  the  well.  We  note  that  the  intensity  and  frequency  of  the 
scarring  contribution  decreases  with  increasing  k^L  as  shown  in  Refs.  [53]  and  [54], 
Thus,  in  the  limit  k0L  — >  oo,  Eq.  (2.31)  applies;  but  the  error  in  (2.31)  depends  on 
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Figure  2.4:  Propagation  directions. 

the  shape  of  the  well.  For  the  time  being  we  will  proceed  on  the  assumption  that  the 
estimate  given  by  Eq.  (2.31)  can  be  used  for  the  denominator  of  (2.24).  Although 
we  will  End  that  (2.31)  works  well  for  the  chaotic  shape  shown  in  Fig.  2.1(b),  we 
will  also  argue  that,  in  other  cases,  (2.31)  may  not  provide  as  good  a  model  for 
fluctuations  of  A E. 

To  invoke  property  (ii)  (i.e.,  isotropy  of  velocity  direction),  we  note  that  the 
terms  in  the  sum  in  (2.26)  represent  wave  propagation  directions  that  make  an  angle, 

„  /  run  \ 

6m  =  arcsm  (  —  J  ,  (2.32) 

with  respect  to  the  x  axis;  see  Fig.  2.4. 

We  imagine  that  these  wave-quantized  angles  represent  a  range  of  the  contin¬ 
uous  classical  angles,  where  the  range  for  9m  is 

A6*m  =  -{0m+  i  -  dm- i),  (2.33) 

and,  for  m  =  M ,  we  replace  9m+  \  by  ir/2,  while,  for  m  =  1,  we  use  9m- 1  =  0. 
Invoking  the  classical  orientation  isotropy  of  particle  velocities,  point  (ii)  above, 


26 


the  correspondence  principle  suggests  that  (c^)  is  proportional  to  A 9m,  (c^)  = 
(constant)  Ad m,  where  the  angle  brackets  denote  an  average  over  our  pseudorandom 
fluctuations.  Using  this  in  (2.31),  we  obtain 

M 

<<&>  =  4«.J)(A0m/  y  A em,).  (2.34) 

m'= 1 

(Note  that  the  sum  over  A dmt  is  (7t/2)  —  (6h/2),  rather  than  ir/2.  We  have  chosen  to 
omit  the  angles  0  <  6  <  6\/2  because  the  normally  incident  wave  corresponding  to 
m  —  0  is  ruled  out  by  our  boundary  condition,  ip0  —  0  at  y  —  0,  L  and  x  =  0.  In  any 
case,  this  choice  makes  only  a  small  difference  for  M  3>  1.)  Since  we  view  the  cm 
as  resulting  from  the  sum  of  many  roughly  independent  classical  ray  contributions, 
the  central  limit  theorem  implies  that  cm  will  be  a  Gaussian  random  variable.  Thus 
we  set 

cm  =  ( c2J1/2U  (2.35) 

where  are  independent,  Gaussian,  zero  mean,  unit  variance  random  variables, 

=  ^m,m' i  (2.36) 

with  5mrni  —  1  if  m  —  m!  and  Smjmi  =  0  if  m! . 

Combining  (2.31),  (3.12),  (2.34),  (2.29),  (2.28)  and  (2.24),  we  arrive  at  our 
main  result, 


where 


A  E  = 


Em  £ 

m=  1  UmS, 


5Em  — 


mam 


Ak2B  sinh  (2cemA) 


(2.37) 


(2.38) 
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is  the  contribution  to  the  splitting  due  to  the  m-th  mode  of  the  barrier, 


dm 


2A  9m 
7T  —  9\ 


(2.39) 


is  the  weight  assigned  to  the  angle  9m  in  the  well,  and  the  Gaussian  random  variables 
Cm  satisfy  (2.36). 

When  M  is  large  (i.e.,  kL  3>  1)  the  number  of  terms  in  the  sums  in  (2.37) 
is  large  and  the  denominator  is  close  to  unity  with  relatively  small  fluctuations. 
Although  the  fluctuations  of  the  denominator  from  unity  are  small  for  large  M,  it 
can  be  necessary  to  include  them,  as  they  significantly  contribute  to  the  evaluation  of 
the  fluctuations  of  A E  from  (A E),  which  are  also  small  for  large  M.  As  before,  the 
angle  brackets,  (...),  represent  an  ensemble  average  over  realizations  of  the  random 
set  {^m}.  The  sliding  average  (...)E,e  hi  Sec.  3.1  is  hypothesized  to  approximate 
(...),  if  (...)  is  approximately  constant  over  the  window  width  e  and  if  many  energy 
levels  are  contained  in  the  window. 

Equation  (2.37)  is  a  statistical  model  for  the  pseudorandom  splittings  A E. 
This  model  can  be  used  to  generate  ensembles  of  values  of  A E  via  the  Monte  Carlo 
procedure  of  generating  and  inserting  random  values  for  the  Gaussian  quantities  Cm, 
from  which  the  statistical  properties  of  A E  can  be  numerically  determined,  as  will 
be  illustrated  by  the  example  given  in  Sec.  2.4. 

For  large  M  ^  kL/n,  both  the  numerator  and  the  denominator  in  Eq.  (2.37) 
will  have  relative  fluctuations  from  their  mean  values  that  are  small.  Recall  that  the 
mean  of  the  denominator  is,  by  construction,  one.  Thus  we  expect  that  replacing 
the  denominator  by  one  (i.e.,  neglecting  denominator  fluctuations)  will  make  little 
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difference  in  the  mean  value  of  A E  obtained  from  Eq.  (2.37).  Fluctuations  of  the 
denominator,  however,  can  have  a  significant  effect  when  looking  at  fluctuations  of 
A E  from  its  mean,  as  we  now  discuss.  Say  the  numerator  has  a  large  upward  (or 
downward)  fluctuation.  This  might  occur  because  £T2n  for  some  m  values  happen  to 
be  significantly  above  (or  below)  their  mean  value  of  one.  If  this  is  so,  then  the 
denominator  will  also  have  a  large  upward  (or  downward)  fluctuation,  and  this  will 
mitigate  the  effect  of  the  numerator  fluctuation  on  A E.  Thus  correlation  of  the 
fluctuations  of  the  numerator  and  the  denominator  reduce  the  overall  fluctuations 
of  A  E. 
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fluctuating  component  is  small  compared  to  unity,  and  we  neglect  it.  (We  emphasize, 
however,  that  inclusion  of  this  fluctuation  can  be  crucial  for  a  calculation  of  the 
variance  of  A E.)  In  the  semiclassical  limit  kL  3>  1 ,  M  1,  and  A 0rn  becomes 
small,  allowing  us  to  approximate  the  summation  over  m  in  (2.40)  by  an  integral. 
Thus  (2.40)  becomes 

(AE)  =  ^KI(E/VB,kBA),  (2.43) 

I(E/Vb,  kB  A)  =  -  r  F(9)d6 ,  (2.44) 

7T  Jo 

and  we  recall  E/VB  =  kl/kB. 

Equation  (2.40)  is  plotted  as  the  red  curve  in  Fig.  2.3.  We  next  illustrate  the 
use  of  (2.37)  by  application  to  the  fluctuation  data  shown  in  Fig.  2.2  (b). 

2.4  Monte  Carlo  Evaluation  of  Energy  Splittings 

In  order  to  quantitatively  compare  our  theory  with  the  numerical  data  for 
energy  level  splittings  in  Fig  2.2,  we  define  the  sliding  average  splitting  and  the 
sliding  average  splitting  variance  as 

\Ecy-Eo\<e 

(A  E)e„,,  =  —  y  A  E„,  (2.45) 

Ne° „ 

\E„—Eo\<e 

4B,a,,e  =  l 5—  E  (AE„-(AE}Eo,,f,  (2.46) 

Ne° .<  „ 

where  NE0,e  is  the  number  of  states  such  that  \E0  —  Ea\  <  e  and  we  choose  e  =  V^o- 
This  quantity  is  plotted  as  a  solid  line  in  Figs.  2.2,  2.3  and  2.5.  To  compare  with  our 

semiclassical  Green  function  approach;  an  example  where  this  approach  has  been  employed  in 
another  context  is  the  paper  by  Ullrno  et  al.  [56] 
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result,  Eq.  (2.37),  we  use  Monte  Carlo  simulations.  At  each  energy  level  Ea  plotted 
in  Fig.  2.2,  we  use  (2.37)  to  generate  10000  splitting  values,  A Eai,  i  —  1,2, ...,  10000. 
Similar  to  Eq.  (2.45)  and  (2.46),  for  each  of  the  10000  set  of  Monte  Carlo  data, 
{A Ea.},  we  also  calculate  sliding  average  splittings, 

\Eai-E0\<e 

(A  E)E„.,,i  =  —  Y.  (2-47> 

Ne°’>  X 

and  sliding  average  splitting  variances, 

\Eai-E0\<e 

<E,E0,.,i  =  Jj—  Y  (AB„,  -  <A£)Eo,.y.  (2.48) 

Ne°*  „ 

At  each  Eq,  we  can  calculate  the  ensemble  average  and  variance  of  the  sliding  av¬ 
erage  splitting  and  the  sliding  average  splitting  variance.  In  Fig.  2.5,  we  compare 
these  Monte  Carlo  results  (plotted  in  red)  with  results  from  numerical  solutions  of 
the  Schrodinger  wave  equation  (plotted  in  black);  we  also  compare  these  results 
to  what  (2.37)  would  predict  if  the  denominator  of  Eq.  (2.37)  were  set  to  unity 
(plotted  in  green).  Figure  2.5(a)  shows  that  no  matter  whether  fluctuations  in  the 
denominator  are  included  or  neglected,  the  sliding  averages  of  the  splittings  for 
both  calculations  fall  within  one  ensemble  standard  deviation  of  each  other,  and 
both  agree  well  with  results  from  numerical  solution  of  the  wave  equation.  In  con¬ 
trast,  we  see  from  Fig.  2.5(b)  that  including  the  fluctuations  in  the  denominator 
reduces  the  splitting  variance.  That  is,  correlations  between  the  denominator  and 
numerator  reduces  the  estimated  eigenfunction  to  eigenfunction  variations  in  the 
splitting  energy. 

To  examine  the  effect  of  the  correlation  between  the  denominator  in  Eq.  (2.37) 
and  the  numerator,  we  numerically  calculate  energy  splittings  for  symmetrical  dou- 
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(a) 


(b) 

Figure  2.5:  Comparison  of  (a)  sliding  average  splittings,  (A E),  and  (b)  sliding 
average  splitting  variances,  (Jae,  versus  E0  for  numerical  data  (black),  Eq.  (2.37) 
(red)  and  Eq.  (2.37)  with  its  denominator  replaced  by  one  (green). 
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V  =  0 


Figure  2.6:  Symmetrical  double  wells  of  area  A  separated  by  a  non  perfectly  coupled 
tunneling  potential  barrier  of  width  2A,  barrier  length  L,  height  Vb  and  wall  length 
w. 

ble  well  that  has  the  same  gross  parameters  (A,  Vb,  A,L)  with  Fig.  2.1  but  longer 
wall  at  x  —  0,  see  Fig.  2.6,  and  make  analogous  figures  to  Fig.  2.5  as  Fig.  2.7.  In 
order  to  explain  the  discrepancy  in  Fig.  2.7(b)  between  the  theory  as  expressed  by 
Eq.  (2.37)  and  onr  data  from  numerical  solution  of  the  Schrodinger  equation,  we  now 
re-examine  our  assumption  that  we  can  use  the  approximation  (2.31)  for  fLW  i^dxdy 
in  (2.37).  In  particular,  as  explained  in  the  discussion  following  Eq.  (2.31),  at  finite 
wavelength,  A^1  fLW  i^dxdy  and  A^1  fA  ipydxdy  may  not  be  perfectly  correlated. 
As  an  alternate  hypothesis,  let  us  now  suppose  that  fluctuations  of  A^1  fA  i/jq dxdy 
are  uncorrelated  with  those  of  A~l  fLW  i/jq dxdy .  If  this  is  the  case,  the  fluctuations  of 
the  denominator  in  Eq.  (2.24)  are  uncorrelated  with  fluctuations  of  the  numerator. 
In  this  situation  we  can  choose  to  constrain  the  denominator  to  be  normalized  to 
one,  fLW  t^d/xd/y  =  1.  In  particular,  this  is  consistent  with  our  previous  definition 
(£m)  =  1  and  conforms  with  the  idea  that  fA  ^ dxdy  averaged  over  many  modes 
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E0 


(a) 


(b) 

Figure  2.7:  (a)Analogous  to  Fig.  2.5(a)  but  for  Fig  2.6.  (b)Analogous  to  Fig.  2.5(b) 
but  for  Fig  2.6. 
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should  respect  the  global  normalization  fLW  i^dxdy  =  1  for  each  mode.  Thus  this 
alternate  hypothesis  yields  (2.37),  but  with  the  denominator  replaced  by  one, 

M 

AE  =  Y^  »mSEmem.  (2.49) 

m=  1 

In  Sec.  2.6.3,  we  provide  analytical  support  for  this  and  show  how  a  transition  from 
applicability  of  (2.37)  to  applicability  of  (2.49)  can  take  place  as  a  geometrical  pa¬ 
rameter  is  varied.  Equations  (2.37)  and  (2.49)  result  from  two  opposite  bases,  perfect 
correlation  for  (2.37)  and  zero  correlation  for  (2.49).  As  previously  discussed  in  the 
text  following  Eq.  (2.37),  correlation  reduces  the  fluctuations.  Hence,  we  expect 
the  fluctuation  level  to  lie  between  the  predictions  from  these  two  extremes,  and 
we  regard  the  green  and  red  variance  curves  in  Figs.  2.5(b)  and  2.7(b)  as  predicted 
upper  and  lower  bounds  for  the  fluctuation  level.  Our  data  for  the  two  chaotic 
shapes,  indeed  conform  to  this  expectation,  with  the  fluctuation  level  for  the  shape 
in  Fig.  2.1(b)  being  close  to  the  lower  bound,  while  that  for  the  shape  in  Fig.  2.6  is 
closer  to  the  upper  bound. 

2.5  Green’s  Function  Analysis  of  Sliding  Average  of  Splittings 

We  have  obtained  an  expression,  Eq.  (2.40)  for  the  average  of  the  splittings, 
(A E),  for  chaotic  cavities.  We  have  also  seen  (Fig.  2.3)  that  this  result  agrees 
numerically  with  our  results  for  an  integrable  cavity  of  rectangular  shape.  Here 
we  demonstrate  in  a  formal  analysis  that  our  result  for  (A E)  applies  for  all  cavity 
shapes  independent  of  whether  they  are  integrable,  chaotic,  or  a  mixture  of  chaotic 
and  integrable  in  different  regions  of  phase  space.  Our  analysis  makes  use  of  previous 
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work  on  perimeter  corrections  to  Weyl’s  equation  for  the  density  of  states  [49]. 

We  begin  with  the  Green’s  function  of  the  unperturbed  left  well  (A  — »  oo  in 
Fig. 2.1)  expanded  in  orthonormal  modes  'ipQ  of  the  left  well. 


Ge(x,x')  =  ^ 

a 


E-Ea 


En 


(2.50) 


where  (V2  +  E)Ge  =  S(x  —  x'),  (V2  +  E o-)^  =  0  with  the  appropriate  boundary 
conditions,  and  x  =  (x,y).  According  to  our  perturbation  theory  (Eq.  (2.23))  the 
splitting  for  unperturbed  mode  a  is 

AEa=  [  {V'o  (x,y)AH[^(x,y)]}x=0-dy,  (2.51) 

Jo 


where  the  operator  AH  =  AH  a  ~  A  H$-  Operating  on  (2.50)  with  AH,  setting 
E  =  E0  —  ie,  where  e  >  0,  we  obtain 


b a  /  { \AH G EQ—ie\x' =x\ x=o~  dy  ^  ^ 

Jo 


:A  En 


(2.52) 


(£0  -  £CT)2  +  e2 

For  e  ^  p_1(E0),  where  p_1(E0)  is  the  average  spacing  between  energy  levels  (p(E0) 
is  the  density  of  states),  yet  small  compared  to  E0,  the  right  hand  side  of  (2.52)  is 
the  product  of  np(Eo)  and  the  Lorentzian  sliding  average  of  AEa.  This  follows  from 


£ 


e/7r  . _  .  A 

(£0-£CT)2  +  e2  “P(  0  =  47t’ 


(2.53) 


where  the  right  hand  equality  is  Weyl’s  formula  for  a  well  of  area  A  and  is  indepen¬ 
dent  of  Eq  for  the  two-dimensional  case  we  are  treating.  Equations  (2.52)  and  (2.53) 
yield  the  following  expression  for  the  sliding  average, 


(^-^cr)  Eq,€  J  Eo—ie)%f=x\  x=ody} 


(2.54) 
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Figure  2.8:  Geometry  for  the  Green’s  function. 


For  the  purpose  of  evaluating  the  right  hand  side  of  (2.54),  we  consider  e  to 
be  large  enough  that  waves  with  energy  E0  —  ie  attenuate  to  negligible  values  in  a 
distance  of  the  order  of  the  well  size.  With  this  stipulation,  we  can  replace  GE0-ie 
in  Eq.  (2.54)  by  G^_ie,  where  G^_ie  is  the  Green’s  function  for  the  case  shown  in 
Fig.  2.8,  with  outgoing  waves  as  x  — *  —  oo. 

Making  use  of  the  delta  function  expansions, 


5{x  r*  x')  —  —  /  exp  [ikx(x  —  x')]dkx , 


s(y-y')  =  ^Y. 


m= 1 


miry'  \  rmny 

sin  I  — - —  I  sm  ' 


(2.55) 

(2.56) 


we  obtain 


Gf  =  —  I  dkx^2 

m=  1 


( nmy'\  (rrmy\elkx<yX  x'^ +  Tme  lk*(x+x')  1 

wL  J  ^  ,  sm  {—)  ™  (—)  E_n  +  (y ),]  )  ■  (2-57) 

Noting  that  for  x  >  0,  the  zero  order  Green’s  function  has  the  form  G^)  =  )T)  Gm 
sin  (■ nrmy/L )  exp  (—amx),  the  reflection  coefficient  Tm  is  determined  from  the  bound¬ 
ary  condition, 

i  Ptn^  I  i  Ptn^  I 

,  (2.58) 


G^}  dx 

E 


Q dx 

x=0~ 


c=0+ 
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which  yields 


T-'  _  ^  2i</>77i  {hx ) 

1  m  c  •> 


ar 


where 


0m  (kx)  =  arctan  (  —  ) , 

& X 


Om=  \lk%+  {^pj2  -Eq. 


(2.59) 

(2.60) 

(2.61) 


Inserting  (2.57)  into  (2.62),  and  making  use  of  our  results  for  Ha,s  in  Eqs.  (2.12) 
and  (2.13),  we  obtain 

Re[  i  +  rm(/cx)] 


4  r+°° 

(A Ea)EQ)e  =  ^  y 1  / 

d  —  OO 


dkT— 


ar 


vr  i  sinh  (2amA)  [fc2  +  (hhl)2  -  E 012  +  e2 


(2.62) 


L  > 

In  writing  (2.62)  we  have  neglected  Im\T m{kx)\  which  will  be  valid  for  Eq  3>  e.  In 
this  same  limit  we  may  also  neglect  the  variation  of  Tm{kx)  and  am  in  the  range, 


e  > 


2  fmn\2 


K+{~)  ~E° 


>  -e. 


Thus  we  set  kx  =  kxo  =  \/ E0  —  ( mn/L )2  and  amo  =  \Jk\  +  (m7r/L)2  —  E0  in  (2.62), 
yielding 

4  J1  ™  Rph  4-  r  (h  h 

(2.63) 


=  4  E  ^mo  i?e[l  +  rm(/c 

xo)] 


A  ^  kx o  sinh  (2am,0  A) 
where  we  have  cut  the  sum  over  m  off  at  M  [defined  by  Eq.  (2.30)]  and  dropped  the 

subscript  E0,  e.  Using  our  result  (2.59)  for  Em,  we  finally  obtain 


4L 


M 


(2 A6 


k 2  a 

rvx0Lxrn0 


(2.64) 


^A  “[  V  71  J  kl  sinh  (2amo  A)  ’ 

where  we  have  used  A 6m  =  n/(kxL),  valid  in  the  limit  m  3>  1.  This  result  is  the 
same  as  our  Eq.  (2.40)  derived  for  the  chaotic  shape,  thus  demonstrating  that  it  is 
independent  of  how  the  well  is  shaped,  as  well  as  whether  the  orbits  are  chaotic, 
integrable,  or  mixed. 
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2.6  Discussion  and  Conclusion 


Defining  fluctuating  weights  wm  by 


w 


m 


Eq.  (2.37)  takes  the  form  of  a  weighted  average, 

M 

AE  ^  wm5Em, 

m=  1 


(2.65) 


(2.66) 


where  S Ern  is  defined  in  (2.38).  For  the  case  of  rectangular  wells,  the  unperturbed 
states  (A  — >  oo)  are 


ipo(x,  V )  =  sin  sin  (kx 


— 


where  m  labels  the  vertical  wavenumber,  ky  =  mn/L. 
perturbation  result  (2.24)  yields 


dm),  (2.67) 

Insertion  of  (2.67)  into  the 


A  E  =  SEm. 


(2.68) 


Thus  A E  in  the  chaotic  case  is  a  weighted  average  over  the  tunneling  rates  for 
the  rectangular  well.  This  self-averaging,  done  by  each  individual  chaotic  mode, 
is  responsible  for  the  reduction  of  the  mode-to-mode  tunneling  fluctuations.  The 
larger  the  number  of  m  values  effectively  taking  part  in  the  averaging,  the  lower  the 
fluctuation  level.  Since  this  number  scales  with  M  «  kL/ n,  we  conclude  that,  as 
shown  in  Sec.  2.6.1,  the  fluctuation  level  for  splittings  in  the  chaotic  case  decreases 
with  increasing  kL, 


(A  E) 


((AE  -  (AE))2)1/2 

(A E) 


(kL) 


-1/2 


(2.69) 
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and  that  the  ratio  of  the  fluctuation  level  for  the  chaotic  case  to  that  for  the  inte¬ 
grate  case  has  this  same  predominant  scaling.  Thus  the  difference  between  fluctu¬ 
ation  levels  of  the  chaotic  and  integrable  cases  becomes  large  with  increasing  energy 
(however,  if  E  is  increased,  Vb  may  also  have  to  be  increased,  in  order  to  keep  E /Vb 
less  than  one). 

Equation  (2.66)  also  provides  a  simple  way  of  understanding  our  observation 
that  the  sliding  averages  for  the  chaotic  and  rectangular  well  cases  are  the  same. 
We  first  recall  that  the  weights  wm  given  by  (2.65)  have  averages  corresponding  to 
an  isotropic  distribution  of  incident  plane  waves  on  the  boundary.  Furthermore, 
according  to  Weyl’s  law  for  two  dimensional  billiards,  the  distribution  of  modes 
in  k- space  is  isotropic  and  uniform.  Thus,  if  the  sliding  average  for  the  rectangle 
includes  a  sufficient  number  of  modes  in  the  averaging  window,  then  it  produces  an 
isotropic  averaging,  just  as  in  the  chaotic  case.  Thus,  as  observed  in  Fig.  2.3  and 
demonstrated  by  the  analysis  of  Sec.  2.5,  we  expect  the  chaotic  and  regular  wells  to 
yield  the  same  sliding  average. 

We  remark  that,  from  the  experimental  point  of  view,  due  to  the  short  wave¬ 
length  necessary  for  observing  semiclassical  effects,  the  symmetry  required  for  real¬ 
izing  splitting  statistics  may  be  stringent.  Another,  much  less  demanding,  situation 
is  that  of  tunneling  from  a  single  well  into  a  region  of  outgoing  quantum  waves, 
as  pictured  in  Fig.  2.9.  In  this  case  the  energy  levels  acquire  an  imaginary  part, 
Ea  =  Ei-R>  —  iE„\  where,  for  E^  <  Vb  and  A  sufficiently  large,  Ei-R>  S>  E„  \  so 
that  a  perturbation  analysis  similar  to  that  in  Secs.  2.2  and  2.3  can  be  applied.  The 
result  is  that  the  statistics  of  the  tunneling  escape  rates  {E^}  are  similar  to  those 
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Figure  2.9:  Tunneling  from  a  single  well  into  an  unconfined  region. 

for  the  tunneling- induced  splittings  {A Ea}.  The  statistical  model  for  { E (anal¬ 
ogous  to  Eq.  (2.37))  is  derived  in  Sec.  2.6.2.  The  result  has  the  same  form  as  that 
given  in  Eqs.  (2.65),  (2.66).  Thus  the  subsequent  discussion,  including  Eq.  (2.69), 
also  applies  for  {Ea'1}.  Hence  the  same  chaos  regularization  of  tunneling  is  expected 
to  apply  for  the  escape  rates  in  situations  like  that  shown  in  Fig.  2.9 

In  conclusion  we  have  presented  a  semiclassical  analysis  of  energy  level  splitting 
of  symmetric,  quantum-dot-type,  double-well  systems,  where  the  wells  are  separated 
by  a  tunneling  barrier.  Our  analysis  quantitatively  explains  the  observed  mean 
splittings  and  their  fluctuations.  The  mean  is  found  to  be  independent  of  the  well 
shape  and  independent  of  whether  the  well  orbits  are  chaotic  or  not.  In  contrast, 
the  fluctuation  statistics  are  vastly  different  when  the  well  orbits  are  integrable,  as 
compared  to  when  they  are  chaotic,  with  the  chaotic  case  yielding  much  reduced 
fluctuations  when  the  lateral  barrier  length  is  large  compared  to  a  wavelength. 
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2.6.1  Variance  of  the  Splittings  for  Large  kL  and  Large  ks A 


Here  we  apply  (2.37)  to  obtain  an  analytical  expression  for  the  variance  of 
(A  Ea}  for  large  kL  and  large  Lb  A.  Because  of  the  approximations  that  we  will  use, 
the  result  will  not  be  applicable  for  explaining  the  numerical  results  in  Fig.  2.2  (b), 
and  this  is  why  we  used  the  Monte  Carlo  procedure  in  Sec.  2.4.  Nevertheless,  this 
calculation  is  instructive;  e.g.,  it  clearly  shows  that  the  splitting  variance  relative 
to  the  mean  (A E)  decreases  as  ( kL )_1  for  increasing  kL  (also  see  the  discussion  in 
Sec.  2.6). 


We  begin  by  using  (2.40)  and  (2.42)  to  re-express  (2.37)  as 

1  -p  OL 


A  E  =  (A  E) 


1  +  /T 


(2.70) 


where 


E".i  - 1) 


a  = 


M 


(2.71) 

(2.72) 


(A  E) 

ft  ^  ^  Lmikm  1)> 

m=  1 

with  (A E)  being  the  expression  given  by  (2.40),  and  (a)  =  (/3)  =  0  by  virtue  of 

<e>  =  i- 

Anticipating  that  a  and  (3  are  small  compared  to  one  (i.e.,  (a2),  (/ 3 2)  C  1),  we 
expand  (2.70)  to  obtain 


A E  -  (A E) 
(A  E) 


=  a  —  f3. 


(2.73) 


Squaring  (2.73)  and  averaging  over  realizations  of  the  Gaussian  random  variables 
{£m}  yields  the  following  expression  for  the  variance  a2. 


<7 


M 


(A  EY 


2E" 

771=  1 


“I  2 


(A  E) 


-  1 


(2.74) 
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where  we  have  used  ((£r2n  —  1)(£^,  —  1))  =  2 6mm>.  For  large  M  (i.e.,  large  kL ),  we 
now  attempt  to  approximate  the  summation  of  the  right  hand  side  of  (2.74)  by  an 
integral  over  9.  For  large  M  and  9m  not  too  close  to  n/2, 


see  Fig.  2.4.  Using  this  in  (2.74)  we  obtain 


-2  _2i 

r6*  2 

\  SE{9)  1 

2  d9 

(A  Ey  J 

o  kL  cos  6 

L  (AE)  \ 

n/2 

(2.75) 


(2.76) 


While  the  upper  limit  of  the  integration  in  (2.76)  might  nominally  be  supposed 
to  be  7t/2,  we  have  instead  replaced  it  by  9 *,  because,  due  to  the  term  l/cos0 
in  the  integrand,  the  integral  diverges  logarithmically  at  9 *  =  tt/2.  This  is  an 
artifact  of  our  approximation  (2.75),  which  is  not  accurate  for  small  cos 9m  (e.g., 
it  predicts  /im  — *  oo  as  9m  — *  vr/2).  Since  the  divergence  is  logarithmic,  the  size 
of  the  contribution  to  the  variance  from  the  vicinity  of  6  near  n/2,  can  be  roughly 
estimated  by  appropriately  cutting  off  the  integral  at  d*  slightly  below  n/2.  Based 
on  the  construction  shown  in  Fig.  2.10  ,  we  choose  as  an  appropriate  cutoff 


(2.77) 


where  7  is  of  order  one.  The  result  will  be  insensitive  to  a  precise  choice  of  7.  The 
contribution  from  6  near  7t/2  is  then  estimated  to  be  of  the  order  of 


(2.78) 


where  we  obtain  this  result  by  approximating  cosd  by  (n/2  —  9)  setting  6  =  n/2  in 
5E(9 ),  and  noting  from  Eq.  (2.41)  that  5E( n/2)  =  0. 
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where  A  =  2n/k.  Thus  we  conclude  from  (2.83)  that,  even  when  L  is  many  wave¬ 
lengths  A,  the  log  contribution,  (2.78),  is  not  significant  and  the  predominant  scaling 
of  a  is  as  in  (2.82), 


(A  E) 


{kL)-1'2. 


(2.84) 


2.6.2  Escape  Rate  from  a  Chaotic  Well  to  an  Open  Region 

In  this  section  we  outline  the  analysis  of  the  situation  illustrated  in  Fig.  2.9. 
Again  taking  x  =  0  to  coincide  with  the  left  face  of  the  barrier,  2A  to  be  the  barrier 
width,  and  L  to  be  the  vertical  dimension  of  the  barrier  boundary,  we  write  il>(x,y) 
in  x  >  2A  and  0  <  x  <  2A  respectively  as 

OO 

^(x,y)  =  ^2  An  sin  ^^jetkx’mX ,  (2.85) 

m=  1 
oo 

if>(x,y)  =  J2(E™e~amX  +  FmeamX)  sin  (^),  (2.86) 

m=  1 

where 

,  /I  rmn\2  „  /m7r\2 

kx,m  =  Je- y—j  for  E  >  >  (2-87) 

,  /  rmir\ 2  Z  „  „  (m  7r\2 

A,m  =  (— )  -  E  for  E  <  J  ,  (2.88) 

=  J 2  +  vb  ~  E,  VB  >  E.  (2.89) 

Applying  the  continuity  of  and  difr/dx  at  x  =  2 A  and  at  x  =  0,  we  obtain  the 
boundary  condition  at  x  =  0-, 

+  H[i/>(0~,y)]  =  0l  (2.90) 

x=0~ 


dx 
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where  H  is  defined  in  a  manner  analogous  to  Eqs.  (2.12)  -  (2.14)  with 


^(rn)  _  (gm  -  ikm)  -  (am  +  ikm)  exp  (— 4amA) 
m  (atm  ~  ikm)  +  ( atm  +  ikm )  exp  (-4amA) ' 


(2.91) 


Proceeding  as  in  Sec.  2.2.3,  perturbation  theory  gives 


E-En=± h 


/qL{V,0A^[V,q]}x=0-^ 
Ilw  ^ odxdv 


(2.92) 


where  E  =  AH  =  H  —  H0,  and  Ho  is  defined  by  (2.18).  Taking  the 

imaginary  part  of  Eq.  (2.92),  we  obtain  an  expression  for  the  tunneling  rate, 


Ed)  _  So{^oIm(AH[^o])}x=o-dy 

JLW^odxdy 


(2.93) 


Assuming  that  exp  (— 4amA)  is  small  we  find  that 


Im(H{m)  -  H^m))  «  -—a2mkx,mex p  (-4amA), 

vb 


(2.94) 


which  yields 


E{i)  =  2L  Em  exp  (— 4amA) 

VB  fTAV  ^dxdy 


(2.95) 


where  (2.95)  is  analogous  to  (2.24).  We  can  now  easily  parallel  the  steps  of  Sec.  2.3 
that  lead  to  Eq.  (2.37).  Indeed,  comparing  (2.95)  and  (2.24),  we  can  obtain  (2.95) 
from  (2.24)  by  making  the  following  replacement  in  (2.24): 


sinh2amA  14 


~Oimkx,m  exp  (  4amA). 


(2.96) 


Using  the  replacement  (2.96)  in  Eq.  (2.37),  we  obtain  the  following  statistical  model 
for  the  tunneling  rates  from  a  chaotic  well  to  an  open  region, 


E =  £  wmE%>, 


(2.97) 
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where 


E(j)  =  8^,m«m.exP(-4amA) 


AVI 


(2.98) 


and  wm  is  as  defined  in  Eq.  (2.65).  Note  that  Eq.  (2.97)  has  the  same  form  as 
Eq.  (2.66). 


2.6.3  Model  for  the  Upper  and  Lower  Bounds  on  the  Splitting  Vari¬ 
ance 

In  this  section  we  consider  a  model  which  is  similar  to  that  in  Fig.  2.6,  but 
with  a  modification  that  will  facilitate  analysis.  This  model  is  shown  in  Fig.  2.11. 
The  main  feature  of  this  model  is  the  addition  of  a  thin  horizontal  hard,  thin  septum 
a  distance  L  from  the  bottom  of  the  center  of  the  billiard.  This  septum  separates 
the  region  near  the  vertical  part  of  the  well  boundary  abutting  the  potential  barrier 
(labeled  Region  2  in  the  figure),  from  that  abutting  the  vertical  hard  wall  well 
boundary  segment  (labeled  Region  1  in  the  figure).  Using  this  model  we  now  present 
an  analysis  supporting  our  claim  that,  as  the  parameter,  L/(L  +  L)  ( L  is  defined  in 
Fig.  2.11),  varies  from  1  to  0,  the  splitting  variance  <r^E,  transitions  from  the  lower 
bound,  Eq.  (2.49)  ,  to  the  upper  bound,  Eq.  (2.37). 

Applying  the  random  plane  wave  hypothesis  to  Region2,  we  model  the  statis¬ 
tics  of  the  spatial  averages  over  Region  2  of  f°r  given  modes  as 

('0o ) 2  =  gmd,  M  =  Int  (  —  J  ,  (2-99) 

m=  1  x  ' 

where  the  overbar  denotes  spatial  average,  Y2m= i  Rm  =  L  are  i-i.d.  Gaussian 
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Hard 


Figure  2.11:  Model  billiard  for  the  analysis  in  Sec.  2.6.3. 


variables, 


(tmtm')  =  (2.100) 

and  a  given  random  realization  of  the  is  hypothesized  to  statistically  model  a 
given  mode.  Doing  the  same  thing  for  Region  1,  we  model  the  statistics  of  spatial 
averages  of  i[}q  over  Region  1  for  given  modes  as 

M)i  =  M  =  Int  (  —  )  ,  (2.101) 

rh=  1  V  / 

with  similar  definition  of  /i„,  and  Averaging  over  many  modes  (such  average  are 
denoted  (•••)), 

m-2)  =  <a,  ((^)i)  =  (?)■  (2.io2) 

Since  we  expect  the  model  averages  of  over  any  region  to  be  the  same,  (£2)  =  (£2), 
and  we  define  (£2)  =  (£2)  =  1. 

We  now  adapt  the  additional  model  hypothesis  of  perfect,  model  by  mode, 
correlation  between  the  average  of  over  the  whole  region  A  and  its  average  over 
the  combined  area  of  Region  1  plus  Region  2.  While  this  may  not  really  apply,  it 
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should  yield  a  valid  qualitative  model  for  the  issue  that  we  wish  to  study.  This  gives 


where 


Area  of  Region2 


(2.103) 


(2.104) 


(L  +  L)  Area  of  Regionsl  +  2 
Letting  r  =  L/(L  +  L)  and  following  our  previous  work,  application  of  Herring’s 

formula  gives 


where 


A  E 


It® 


(A  E)  1  +  r/3  +  (1  —  r)y’ 


M 

E 


SE„ 


_  Jm  ( 1  \ 

a  —  ft m  m  ■*■)) 


m- 

M 


ft  —  ^  ^  ftm{£,m  1); 

m=l 

M 

7  —  ^  ftrhdrh  ~ 


(2.105) 


m=l 


(a)  =  (/3)  =  (7)  =  0. 


(2.106) 


Expanding  for  small  a  P  ~  7  <  1,  we  obtain  the  following  expression  for 
the  normalized  splitting  fluctuation  5, 


r  _  A E-  (A E)  „ 

s=  m  = 

The  normalized  splitting  variance  is  thus 


r(3  +  (1  —  r)y. 


(2.107) 


(<52>  9*  («2)  -  2 r(a(3)r2{(32)  +  (1  -  r)2(y2),  (2.108) 

where  we  have  used  (ay)  =  (P'j)  =  0,  reflecting  the  assumption  that  (£m£m)  —  0  for 
all  m  and  m. 
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Say  we  keep  k,  L  and  A  fixed  and  vary  L,  which  is  equivalent  to  keeping  a 
and  (3  fixed  and  varying  r.  How  does  (52)  change  as  r  varies? 

From  Eq.  (2.78)  of  Sec.  2.6.1 


(n 

Thus  for  kL  ~  kL  S>  1,  we  have 
expression  for  (y2),  gives 


L  ln(/cL/27ry) 
L  ln(/cL/27T7) 


(2.109) 


(y2)  =  ( L/L)(/3 2),  which,  when  used  in  our 


(52)  =  r{(a  —  /3)2)  +  (1  —  r)(a2),  1  >  r  >  0.  (2.110) 


Thus  (S2)  varies  linearly  from  its  largest  value,  (a2)  at  r  —  0  (corresponding  to 
Eq.  (2.49)),  to  its  smallest  value,  ((a  —  f3)2)  at  r  —  1  (corresponding  to  Eq.  (2.37)). 
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Chapter  3 


Statistical  Model  of  Short  Wavelength  Transport  Through  Cavities 
with  Coexisting  Chaotic  and  Regular  Ray  Trajectories 
3.1  Introduction 

In  principle,  for  a  given  configuration,  properties  of  wave  systems  are  com¬ 
pletely  determined,  and  thus  are  not  random.  However,  at  short  wavelength,  these 
properties  can  be  very  sensitively  dependent  on  small  configurational  changes  or 
changes  of  the  free  space  wavelength.  If  the  configuration  or  free  space  wavelength 
is  regarded  as  slightly  uncertain  within  some  small  range  and  the  wave  properties 
vary  wildly  in  this  range,  then  a  statistical  approach  may  be  warranted.  This  type 
of  approach  was  originally  introduced  by  E.  Wigner  in  reference  to  the  energy  levels 
of  large  nuclei  [6,  7,  8,  9,  1,  2],  and  later  employed  to  study  classically  chaotic  quan¬ 
tum  systems  [1,  3].  Here  we  focus  on  quasi- two-dimensional  microwave  cavities  and 
quantum  dots  which  couple  to  an  external  environment  through  suitable  openings 
(called  ‘leads’  or  ‘ports’).  The  statistical  properties  in  chaotic  cavities  with  external 
connections  have  been  well  studied  using  various  approaches,  e.g.,  the  ‘Poisson  Ker¬ 
nel’  [18,  19,  20,  21,  14]  or  the  ‘Random  Coupling  Model’  (RCM)  [57,  58,  59].  The 
RCM  (employed  in  this  chapter)  focuses  on  impedance  matrices  (related  to  scatter¬ 
ing  matrices  through  an  elementary  transformation)  and  replaces  the  eigenfunctions 
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and  eigenenvalues  in  the  impedance  formula  by  suitably  choosen  random  quantities. 
Past  work  has  shown  that  the  RCM,  and  , equivalently,  the  Poisson  Kernel  yield 
results  that  agree  well  with  statistical  data  obtained  from  experiments  and  numer¬ 
ical  computations  on  microwave  cavities  [12,  13,  14,  15].  However,  in  general,  such 
systems  may  have  not  only  either  all  chaotic  or  all  regular  orbits,  but  also  typically 
have  a  mixture  of  coexisting  chaotic  and  regular  orbits.  We  called  such  systems 
‘mixed’.  The  statistical  properties  of  impedance  matrices  in  mixed  systems  is  the 
subject  of  this  chapter. 

For  specificity  we  focus  on  a  particular  mixed  system,  a  ‘mushroom’  cavity 
(Fig.  3.1(a))  [60],  which  has  a  clearly  divided  phase  space  [61] . 1  For  most  modes 
of  this  system,  we  find  that  it  is  possible  to  separate  them  into  two  classes,  regular 
and  chaotic  (this  may  not  hold  for  other  systems).  Using  this  separation,  we  de¬ 
compose  the  impedance  formula  into  chaotic  and  regular  parts.  We  then  derive  the 
probability  distribution  associated  with  the  chaotic  part  of  the  impedance,  while, 
for  the  regular  part  we  utilize  exact  (numerically  calculated)  or  approximate  the¬ 
oretical  eigenmodes.  To  test  our  theory,  we  numerically  solve  for  eigenvalues  and 

eigenfunctions  of  our  mushroom  cavity  and  insert  them  into  the  exact  formula. 

1More  generic  systems  can  display  infinite  hierarchies  of  KAM  island  chains  encircling  other 
KAM  island  chains  with  chaos  intermixed.  This  type  of  intermingling  of  chaotic  and  nonchaotic 
orbits  on  all  scales  is  not  present  in  the  mushroom  billiard  where  non-smooth  shape  is  designed  to 
yield  a  clear  division  between  chaotic  and  regular  phase  space  region.  Our  motivation  in  using  the 
mushroom  shape  is  that  the  simplicity  provided  by  its  clear  division  of  chaotic  and  regular  phase 
space  allows  a  potentially  simpler  theory.  We  hope  that  our  work  can  serve  as  a  basis  for  future 
study  applicable  to  the  case  of  generic  phase  space  structure 
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This  chapter  is  organized  as  follows.  In  Sec.  3.2  we  review  the  impedance  for¬ 
mula  in  two  dimensional  cavities,  introduce  the  random  coupling  model,  generalize 
the  RCM  to  mixed  systems,  introduce  the  mushroom  cavity  (an  example  of  a  mixed 
system),  and  apply  onr  generalized  RCM  to  this  cavity.  In  Sec.  3.3  we  numerically 
calculate  the  impedance  matrix  of  the  mushroom  cavity  and  compare  the  numeri¬ 
cal  results  with  results  from  our  statistical  theory.  Conclusions  and  discussion  are 
presented  in  Sec.  3.4. 

The  general  problem  of  wave  properties  of  systems  whose  ray  equations  have 
a  mixed  phase  space  was  first  addressed  by  Berry  and  Robnik  [43]  who  studied 
the  spectra  of  mixed  closed  systems.  Subsequently,  many  other  researchers  have 
investigated  spectra  and  wavefunctions  of  closed  systems  with  mixed  ray  orbit  phase 
space  (e.g.,  [27,  64,  65]).  The  problem  of  characterizing  the  input/output  properties 
of  mixed  open  systems,  however,  has,  to  our  knowledge,  been  addressed  relatively 
little  [66,  67,  68], 

3.2  Review  of  Theory 
3.2.1  Impedance  of  a  cavity 

In  the  presentation  that  follows,  we  consider  the  context  of  electromagnetic 
waves.  However,  we  emphasized  that,  with  appropriate  notational  changes,  these 
considerations  apply  equally  well  to  quantum  waves,  acoustic  waves,  elastic  waves, 
etc. 

We  consider  a  vacuum-filled,  quasi-two-dimensional  (vertically  thin)  microwave 
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cavity  with  cavity  height  h  and  M  ports  as  shown  in  Fig.  3.1.  We  denote  the  two 
dimensional  interior  of  the  cavity  by  fl  G  M2.  If  the  frequency  is  not  too  high  (i.e., 
the  wavelength  is  greater  than  2 h),  then  only  vertical  electric  fields  are  excited  inside 
the  cavity, 

E  =  Ez(x,t)z,  (3.1) 

where  x  G  11  is  a  two  dimensional  position  vector.  The  surface  charge  density  on 
the  bottom  plate  of  such  a  cavity  is  ps  =  —e0Ez,  and  the  voltage  difference  between 
the  two  plates  is 

Vr(x,t)  =  hEz(x,t).  (3.2) 

The  surface  current  density  on  the  bottom  plate  is  related  to  the  magnetic  field,  H, 
which  is  perpendicular  to  E,  by 

Js  =  H  x  z.  (3-3) 

We  assume  that  the  fields  are  excited  by  M  localized  current  sources,  which  inject 
surface  charge  density  on  the  bottom  plate 

M 

Ps(x,t)  =  Y  tjitfujix),  (3.4) 

3= 1 

where  u3(x)  is  the  normalized  profile  function  of  port  j,  f  d2xuj(x )  =  1,  and  we 
regard  Eq.  (3.4)  as  modeling  the  currents  induced  by  the  transmission  line  fed  ports 
shown  in  Fig.  3.1.  With  Eq.  (3.3),  the  continuity  equation  for  the  surface  charge 
can  be  written  as 

q  M 

-(-e0Ez)  +  \7  ■  (H  x  z)  =  ps  =  Y  1iur  (3-5) 

i=i 


54 


(a) 


Coaxial 


Ez 


(b) 

Figure  3.1:  (a)  Top  view  of  the  quasi-two  dimensional  cavity  coupling  with  M  —  2 
ports  (fed  by  coaxial  transmission  lines),  where  the  region  interior  to  the  cavity 
is  denoted  fl  (b)  Side  view  of  the  cavity  at  a  port.  In  some  previous  works,  a 
mushroom  billiard  similar  to  that  in  (a)  was  used  [60],  but  the  billiard  section  below 
the  quarter  circular  cap  being  a  rectangle  of  width  pQ.  This,  however,  introduced 
neutrally  stable  ray  orbits  that  bounce  back  and  forth  horizontally  between  the 
vertical  walls  of  the  rectangle.  By  using  the  above  triangular  bottom  part  (as  in 
Ref.  [63])  (a)  we  avoid  the  non-generic  effects  of  such  orbits. 
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Differentiating  Eq.  (3.5),  using  Faraday’s  law,  Vx£  =  ~n0dH/dt,  and  expressing 

Ez  by  Eq.  (3.2),  we  obtain 

1  02  M  d 

Vt  ~  WT  =  hy „  ^  u^Ij,  (3.6) 

3= 1 

where  c  =  1/y^oeo  is  the  speed  of  light.  Assuming  that  Vr(x,t )  =  Vr(x)ejwt, 
Jj(t)  =  Eq.  (3.6)  can  be  rewritten  as 

M 

(V2  +  fc2)ET  =  trjkhrto^Ujij,  (3.7) 

j=i 

where  k  =  lu/c,  and  r/o  =  y7 /io/eo  is  the  free  space  impedance. 

We  expand  Vt  in  the  basis  of  the  eigenfunctions  of  the  closed  cavity,  i.e., 

OO 

VT  =  ^cn0n,  (3.8) 

n=l 

where  (f)n  satisfies  the  Helmholtz  equation  with  Dirichlct  boundary  condition  and  a 
proper  normalization  condition,  i.e., 


(v2  +  k2n)<t>n(x)  =  0  fen, 


(3.9) 


4>n{x)  =0  X  G  d n, 


kfyjd  %  &ij i 


(3.10) 

(3.11) 


and  we  order  the  mode  labeling  according  to  the  convention,  k2l+ ,  >  k 2.  Inserting 
Eq.  (3.8)  into  Eq.  (3.7),  multiplying  0m(x)  and  integrating  over  D,  we  obtain 


M 


-jkhrj  o  ^ 


l=i 


(' uj(t)m)ij 

k2  —  k2  ' 


(3.12) 


where  (•••)  =  fQ  ■  ■  ■  d2x.  The  voltage  at  port  ?’  is  defined 


as 


Vi  =  («Wt>, 


(3.13) 
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where  the  port  voltages  R  are  expressed  in  phaser  form,  Vt  =  Vle3Ujt.  Using 
Eqs.  (3.8),  (3.12)  and  (3.13),  we  obtain 


M 

3= 1 

where  the  i ,  j  element  of  the  impedance  matrix  Z  is  given  by 

CMn)  ( UjCpn ) 


(3.14) 


=  -jkhr]0  X 

n=l 


k2  -  kl 


(3.15) 


Equation  (3.15)  states  that,  in  a  lossless  cavity,  the  impedance  is  purely  imagi¬ 
nary,  since  the  eigenfunctions  for  Eqs.  (3.9)  and  (3.10)  are  real.  It  also  states  that,  if 
we  know  all  the  eigenfunctions  and  eigenvalues  of  the  closed  cavity,  we  can  calculate 
the  matrix  elements  of  Z  exactly.  Note  that  ( U{(j)n )  — »  0  as  the  port  size  becomes 
much  greater  than  several  wavelengths.  Thus,  the  infinite  sum  in  Eq.  (3.15)  can  be 
replaced  by  a  finite  sum,  i.e., 


N 


Zij  =  -jkhrjo  ^ 


n= 1 


k2-k2n  : 


(3.16) 


where  N  satisfies  the  condition,  2 ti /k^  <C  (size  of  ports).  For  systems  that  are  large 
compared  to  a  wavelength  (27 r/k)  and  may  have  some  uncertainty  in  their  specifi¬ 
cation,  it  is  often  of  practical  interest  to  dispense  with  the  necessity  of  numerically 
calculating  all  N  eigenfunctions  and  to  instead  look  for  a  statistical  description.  The 
later  will  be  our  goal. 


3.2.2  Random  Coupling  Model 

The  Random  Coupling  Model  (RCM)  treats  the  case  where  typical  ray  orbits 
are  all  chaotic  and  is  based  on  the  supposition  that,  in  the  short  wavelength  limit, 
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the  statistical  properties  of  the  impedance  of  a  chaotic  cavity  can  be  obtained  from 
Eq.  (3.16)  by  replacing  k'^  and  ( Ui<f>n )  by  suitable  random  variables. 

According  to  the  Weyl’s  formula  [39]  for  a  two  dimensional  cavity  of  area  A, 
the  mean  spacing  between  two  adjacent  eigenvalues,  —  k^_ i,  is  An /A,  i.e., 

A  =  (kl  -  kl _,)  =  E.  (3.17) 

References  [6,  7,  8,  9,  1,  5]  state  that  the  normalized  eigenvalue  spacing,  sn  = 
(k^  —  fe^_1)/A,  of  a  time-reversible  chaotic  system  has  similar  statistical  properties  to 
the  spacings  of  the  eigenvalues  of  large  matrices  randomly  drawn  from  the  Gaussian 
Orthogonal  Ensemble  (GOE)  of  random  matrices  with  unit  mean  eigenvalue  spacing. 
In  this  chapter,  our  eigenfunctions  are  always  real,  as  appropriate  to  time  reversible 
systems,  and,  henceforth,  GOE  is  automatically  assumed  when  we  mention  random 
matrices. 

Berry  [44]  argues  that  the  wavefunction  at  any  point  in  a  chaotic  billiard  has 
similar  statistical  properties  to  a  random  superposition  of  many  plane  waves, 

4>n(x)  ~  Re  exP  •  x  +  i/3j)  j>  ,  J>  1,  (3.18) 

where  it  is  assumed  that  x  is  not  too  close  to  the  billiard  boundary,  the  wavenumber 
kn  is  fixed,  but  propagation  directions  e},  amplitudes  aq,  and  phases  / 3j  are  random 
variables.  To  be  more  specific,  directions  and  phases  are  uniformly  distributed  in 
[0,  2:/r],  and  all  amplitudes  have  the  same  distribution.  By  the  central  limit  theorem, 
for  J  3>  1,  4>n(x)  evaluated  at  the  point  a:  is  a  Gaussian  random  variable  with  zero 


mean,  and  its  variance  can  be  determined  by  the  normalization  condition,  i.e., 


(3.19) 


which  implies 

E{^}  =  (3.20) 

The  probability  distribution  function  of  the  overlap  integral  ( tq0n )  is  Gaussian 
with  expectation  value  zero  (since  <i>n  is  a  Gaussian  with  expected  value  zero),  and 
by  Eq.  (3.18)  the  variance  of  ( Ui4>n )  is 

i  r 2n  in 

E  {Mn)2}=jJo  u{kn)\\  (3.21) 

where  kn  =  (kn  cos  6,  kn  sind),  and  u[kn )  is  the  Fourier  transform  of  the  prohle 
function  u(x), 

u(kn)  —  j  d2xu(x)  exp  (— ikn  ■  x).  (3.22) 

Note  that,  the  variance  of  ( Ui<pn )  depends  on  the  eigenvalue  k2n  through  Eq.  (3.22) 
where  \kn\  =  kn.  If  2n/kn  3>  (size  of  the  port),  the  prohle  function  of  the  port  can 
be  approximated  by  a  delta  function,  i.e.,  ( Ui<f)n )  =  0n(Tj);  if  2tt /kn  is  comparable 
to  the  port  size,  we  need  to  consider  the  variations  of  <f>n  over  the  ports.  Eventually, 
for  short  enough  wavelength  we  have  E {(u.j0n)}  — *  0  as  kn  — >  oo. 

For  an  M  port  system,  we  need  to  consider  the  same  wavefunction  at  different 
positions;  e.g.,  if  27 r/k  3>  (size  of  the  port),  for  two  ports  located  at  ay  and  x3 ,  we 
need  to  consider  («j0n)  =  0n(^i)  and  ( Uj(j)n )  =  4>n{xj ),  which  are  not,  in  general, 
independent,  although  independence  can  be  approximately  assumed  if  the  ports  are 
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many  wavelengths  apart.  In  the  RCM,  we  build  in  this  relation  by  writing 


$n  =  [(uM,  ■  ■  • ,  {umK)}'1  = 


VA 


W  r 


(3.23) 


where  the  \j\[A  factor  is  based  on  the  expectation  value  of  and  w n(n  = 
1,2, ... ,  N )  is  an  M-dimensional,  zero  mean,  standard  Gaussian  random  vectors 
whose  covariance  matrix  may  have  nonzero  non-diagonal  elements  reflecting  corre¬ 


lation  between  nearby  ports.  We  can  rewrite  the  impedance  matrix  as 

A  N 

2,=  -jkhmj 


A  N  rp 

A  x  - 

'47T  '  k2  —  k2 

n= 1  n 


(3.24) 


where  we  have  used  Eq.  (3.17)  to  replace  A. 

In  the  case  of  identical  transmission  line  inputs  that  are  far  enough  apart,  we 
can  neglect  correlations  between  the  ports  and  the  covariance  matrix  of  w„  is  Imxm] 
i.e.,  E{wlw J')  =  Sij  for  i ,  j  =  1,2,...,  M.  In  this  case,  we  introduce  the  normalized 


reactance  matrix, 


7 T  '  ^ 


(3.25) 


71  n  &  ~  K' 

where  k 2  =  k2 / A  and  the  mean  spacing,  k2  —  k2_x,  between  normalized  eigenvalues 
is  one.  In  this  case  the  impedance  matrix  becomes 


Z  =  j 


.  kh  r/o  „ 


(3.26) 


Note  that  the  normalized  reactance  matrix,  H,  is  independent  of  all  system  specihc 
information,  such  as  the  cavity  shape,  area,  etc.;  namely,  it  is  universal  for  all  chaotic 
cavities  with  widely  separated  ports. 
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3.2.3  Impedance  in  Mixed  Systems 


For  a  generic  two  dimensional  billiard,  both  regular  and  chaotic  phase  space 
regions  coexist,  and  we  call  such  a  system  mixed.  Percival’s  conjecture  [42]  states 
that  semiclassical  eigenmodes  in  mixed  systems  live  either  in  regular  or  chaotic  re¬ 
gions.  Our  numerical  computations  support  this  conjecture  (see  Figs.  3.2  and  3.3). 
At  short  wavelength,  the  number  of  regular  and  chaotic  eigenstates  can  be  approx¬ 
imately  counted  by  the  Partial  Weyl  law  [69], 


Nr(k2)  =  ^k2  +  0(k),  (3.27) 

47T 

where  V  =  R  denotes  regular  trajectories  and  V  =  C  denotes  chaotic  trajectories, 
Ar/A  is  the  ratio  of  the  phase  space  volume  occupied  by  T,  and  Ar  is  given  by 


At  =  [  d2x—  [  d9(r(x,9).  (3.28) 

Jn  ^  Jo 

Here,  Cr {x,9)  is  the  characteristic  function  of  V  at  (x,9),  i.e.,  Cr (x,9)  =  1  if  the 
trajectory  running  through  x  at  9  angle  belongs  to  T  and  Cr(x,  9)  =  0,  otherwise. 

Following  the  above  approach,  we  decompose  (3.16)  into  the  contributions  Zr 
and  Z c  to  the  impedance  from  the  regular  eigenmodes  and  chaotic  eigenmodes,  as 
follows, 

Z  =  Zr  +  Z  c,  (3.29a) 


and 


Nr 

ZR,ij  =  —jkhr]0 


NC 

Zc,a  =  -jkhr)  o  ^2 

C 


(UjCpr) 
Jc2  —  tc 2 

rv  rvr 

(Ui(f>c) 

k2  -  k2c 


(3.29b) 

(3.29c) 
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where  0r(0c)  denotes  regular  (chaotic)  wavefunctions,  r  —  1,2, ,  NR(c  =  1,2, ... ,  Nc), 
and  Nr  +  Nc  =  N. 

The  semiclassical  wavefunction  distribution  for  chaotic  eigenfunctions  in  mixed 
systems  can  be  described  by  the  so-called  Restricted  Random  Wave  Model  [70], 


PiM=^mexp 


2cr2  (d?)  J  : 


(3.30) 


where 


a2(f)  =  ^4):  J  d0(c(x,9). 


(3.31) 


In  a  two  dimensional  pure  chaotic  cavity,  a2  =  1/A  is  independent  of  x. 


The  statistics  of  k2  in  mixed  systems  is  hypothesized  to  be  similar  to  the 
statistics  of  k2  in  chaotic  systems,  but  the  mean  of  the  spacing  between  chaotic 
eigenvalues,  k2+l  —  k2,  is  given  by  4 n /Ac,  as  opposed  to  An /A  in  the  purely  chaotic 
case.  Thus,  the  statistics  of  the  chaotic  normalized  reactance  in  mixed  systems 
should  be  identical  to  the  statistics  of  the  normalized  reactance  in  chaotic  systems. 

We  do  not  expect  to  find  explicit  universal  statistics  for  the  regular  eigenfunc¬ 
tions  (f)r  as  they  are  dependent  on  the  cavity  shape.  However,  the  regular  normalized 


reactance  in  mixed  systems  is  always  Lorentzian  distributed  (see  Appendix  E). 


3.2.4  Mushroom  Billiard 

The  mushroom  billiard  [60,  71]  was  first  introduced  by  Bunimovich.  Since  the 
cap  of  the  mushroom  is  a  quarter  circle,  there  are  orbits  that  never  leave  the  cap 
region  and  are  the  same  as  the  orbits  in  a  complete  quarter  circle  billiard  having 
the  same  radius  R  (see  Fig.  3.2(a)).  These  orbits  are  tangent  to  a  circular  caustic 
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with  a  radius  CJr.  If  the  caustic  radius  Cr  >  p0,  (see  Fig.  3.1)  this  orbit  is  trapped  in 
the  cap,  and  is  integrable.  There  are  also  chaotic  orbits  that  travel  throughout  the 
whole  billiard  (Fig.  3.2(b)),  visiting  both  the  cap  region  and  the  triangular  region 
below  the  cap.  Thus,  the  mushroom  billiard  is  an  example  of  a  mixed  system. 

The  eigenmodes  of  the  Helmholtz  equation  in  a  quarter  circle  with  radius 
R  can  be  described  by  two  quantum  number,  (■ m,n )  -B-  r,  and  the  corresponding 
eigenfunction  is 

4>r  =  4>mn(Pi  #)  =  -A fmnJm(kmnP)  Sin  Uld ,  (3.32) 


with  normalization  constant 


2V2 


(3.33) 


y/nRJ'JkmnR)  ’ 

and  (j)mn  =  0  outside  the  quarter  circle.  Here  Jm  is  m-th  order  Bessel  function  of 
the  first  kind,  krnn  is  the  eigenwavenumber  such  that  kmnR  is  the  n-th  zero  of  Jm. 

To  relate  the  quantum  eigenmodes  to  the  classical  motion  [72],  we  first  define 
the  classical  probability  distribution  for  position  p, 


Pgl(p)  v^r2  -c;Vp2  -cf 

where  PcL(p)dp  represents  the  fraction  of  time  a  classical  trajectory  spends  in  the 
interval  dp  at  p,  R  >p>  Cr.  The  classical  caustic  radius  Cr  is  defined  in  terms  of 
the  angle  of  incidence  (j)  that  the  trajectory  makes  with  respect  to  the  boundary  at 
R ,  Cr/R  =  sin  (j).  The  analogous  caustic  radius  Cr  from  the  wavefunction  (3.32)  is 
identified  by  equating  the  Bessel  function  order  to  its  argument, 


Cr 


Rr 


m 


-R. 


(3.35) 
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(a) 


(b) 


Figure  3.2:  (a)  Two  regular  orbits  with  slightly  different  initial  conditions,  (b)  Two 
chaotic  orbits  with  slightly  different  initial  conditions. 
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Figure  3.3:  (a)  Magnitude  squared  of  the  n  «  10,  002-th  eigenmode  (regular)  and 
kn  253.496413.  (b)  n  «  10,  003-th  eigenmode  (chaotic)  and  kn  f=a  253.501722. 


For  eigenmodes  with  Rmn  <  Po,  the  classical  orbit  in  the  full,  quarter-circle  billiard 
will  travel  to  the  root  of  the  mushroom  so  the  orbit  in  the  mushroom  is  no  longer 
integrable,  and  the  corresponding  modes  in  (3.32)  are  not  present  in  our  system. 
Thus,  we  can  approximate  (3.29b)  using  the  quarter  circle  eigenfunctions  (f)mn  given 
by  Eq.  (3.32), 


Zij,R  =  -jkh'Po 


£ 

m,n 


k2  -  kl 


(3.36) 


Po<R  mn  <R 

In  order  to  apply  the  RCM  for  the  chaotic  contribution  to  the  mushroom 
cavity,  we  need  the  statistics  of  k2  (the  eigenvalues  of  the  chaotic  modes)  and  4>c(x) 
(the  corresponding  eigenmodes).  The  distribution  of  k2  is  taken  to  be  the  same  as 
that  of  the  eigenvalues  of  a  random  matrix  with  same  mean  spacing  Ac  =  (k2+]  — 
kl)  =  At: /Ac-  Using  Eq.  (3.28),  we  can  calculate  the  equivalent  chaotic  area  of  the 
mushroom  cavity, 


Ac  = 


v/3 


2  Po  +  2 


PoyJ R?  -  pi  +  R2  arcsin 


(3.37) 


To  develop  a  random  coupling  model  in  a  mixed  system,  we  need  to  rewrite 
Eq.  (3.23)  as 


<Fn  =  Qwn, 


(3.38) 


where  Q  is  a  M  x  M  diagonal  matrix,  which  describes  the  classical  chaotic  proba¬ 
bility  at  each  port 

Q\  =  /  ui{x)o1{x)(i2x^  (3.39) 

Jn 

where  &(x)  has  been  dehned  in  Eq.  (3.31).  Thus  in  the  case  where  all  transmission 
lines  are  identical,  the  chaotic  contribution  to  the  impedance  matrix  (3.29c)  can  be 
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written 


Zc,ij(k2)  =  J^AcQuQnEtr  (3.40) 

Figures  3.3(a)  and  3.3(b),  respectively,  show  representative,  numerically  com¬ 
puted,  regular  and  chaotic  eigenfunctions.  These  figures  and  others  (not  shown) 
demonstrate  that,  consistent  with  Percival’s  conjecture  [42],  the  eigenfunctions  con¬ 
centrate  either  in  the  regular  or  chaotic  phase  space  regions  thus  justifying  the 
decomposition  (3.29).  We  next  test  the  statistics  predicted  by  Eq.  (3.40)  by  com¬ 
parison  with  direct  numerical  computations  on  our  mushroom  billiard  example. 

3.3  Numerical  Experiment 

In  order  to  test  our  theory  for  the  impedance  in  mixed  system,  we  numerically 
solve  the  Helmholtz  equation  for  its  eigenfunctions  and  eigenenvalues  to  calculate 
Eq.  (3.29)  and  compare  with  our  statistical  model,  Eqs.  (3.36)  and.  (3.40).  We  use 
about  10,000  eigenmodes  for  the  sum  in  Eq.  (3.16).  For  our  numerical  eigenmode 
solutions,  we  use  the  scaling  method  introduced  by  Vergini  and  Saraceno  [73,  63] 
which  facilitates  relatively  fast  solutions.  More  detail  of  this  numerical  technique 
is  described  in  Appendix  B.  It  has  already  been  shown  that  this  method  yields 
accurate  results  for  the  eigenmodes  of  the  Mushroom  billiard  [71].  We  use  a  =  3/4 
(see  Fig.  3.1(a))  rather  than  the  value  a  =  2/3  employed  in  Ref.  [71],  in  order 
to  allow  application  of  Steed’s  Method  [74]  for  efficient  evaluation  of  the  Besssel 
function. 

After  solving  for  all  eigenmodes,  we  classify  these  eigenfunctions  by  examining 
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the  magnitude  of  their  normal  derivative  as  a  function  of  the  boundary  coordinate  s 
(see  Fig.  3.4).  By  this  means  we  can  associate  all  our  numerically  calculated  regular 


(b) 


Figure  3.4:  (a)Regular  eigenmode,  0m, 3(0?),  in  (b) Corresponding  magnitude  of 
the  normal  derivative  of  (ftu,3(x)  versus  s. 

eigenmodes  with  one  of  the  analytically  predicted  approximate  eigenmodes  (3.32). 
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Moreover,  we  have  also  compared  the  regular  eigenfunctions  and  eigenenvalnes  de¬ 
termined  by  onr  numerical  solutions  with  the  approximate  analytic  solutions;  they 
agree  well.  Thus,  the  regular  contribution  to  the  impedance  matrix  (3.29b)  is  very 
well-approximated  by  Eq.  (3.36)  with  onr  approximate  analytic  regular  eigenfunc¬ 
tions  (3.32).  [Alternatively,  one  can  also  characterize  the  regular  contribution  to 
Z  in  a  more  universal  manner,  independent  of  specific  geometry,  as  described  in 
Appendix  E.] 

Onr  first  goal  is  to  test  onr  statistical  model  for  the  chaotic  contribution  to 
Zc  =  Z  — Z  R,  where  onr  model  requires  only  simple  system  information  (cavity  area, 
phase  space  distribution)  rather  than  all  numerical  eigenfunctions.  For  simplicity, 
we  choose  all  ports  to  be  identical,  uncorrelated  and  point- like,  i.e.,  Ui(x)  =  5(x—Xi ); 
thus,  Qu  =  (j(xi)  and  Eq.  (3.16)  becomes 

Zij{k 2)  =  jkhrio&jik2),  (3.41) 


where 


N 


C  (U2\  _  ^nix^C^niXj 

2 L  k2-ki 


n=l 


(3.42) 


and  we  similiarly  define  and  £r- 

We  choose  the  cutoff  Nc  =  N  x  Ac/A  =  2k2/ Ac.  With  this  definition,  the 
expectation  value  of 


2k2  /  Ac 


ic,ij[k2)  =  y 


4>c{Xi)(j)c{Xj 


(3.43) 


c=  1 


k2  —  k/  ’ 

is  zero  since  we  expect  equal  number  of  k2  such  that  k2  >  k2  and  k2  <  k2.  Our  goal 
is  to  find  the  probability  density  functions  of  £c,ij  if  we  randomly  choose  a  k2  (see 
Fig.  3.5). 
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Figure  3.5:  Numerical  calculation  of  £c,u  (red  triangle)  and  £c,ij  (black  square)  in 
the  mushroom  cavity  vs.  energy  ( k 2). 

We  use  a  Monte  Carlo  method  to  generate  realizations  of  Eq.  (3.43).  In  each 
realization,  we  generate  k\ ,  k%, . . . ,  k%c  by  calculating  the  eigenvalues  of  a  GOE 
random  matrix  and  unfold  the  spectra  [1]  such  that  the  mean  spacing  is  47t /Ac', 
we  also  generate  (fa (£*),  fa (£,)),  (fa(xi),  fa(xj)), . . . ,  {(pNc(xi),  <j>Nc(xj))  according 
to  Eqs.  (3.30)  and  (3.31);  then,  we  calculate  £c,ij  at  each  value  of  fc2;  finally,  we 
construct  a  probability  density  function  for  £c,ij-  After  Nr  realizations,  we  have  Nr 
probability  density  functions  for  £c,ij >  be,  pn{£),n  =  1,  ■  ■  ■  ,Nr.  We  then  calculate 
the  mean  and  variance  of  the  probability  density  at  each  £,  i.e, 


(3.44) 


(3.45) 


We  also  calculate  Eq.  (3.43)  numerically  for  different  port  positions  from 


the  numerically  determined  eigenfunctions  and  eigenvalues  and  compare  with  our 
statistical  model  Monte  Carlo  method  (see  Fig.  3.6).  Our  statistical  model  of 
impedance  in  different  port  positions  is  the  statistical  model  of  the  same  normal- 


70 


ized  impedance  (Eq.  (3.25))  with  a  position  dependent  factor,  AcQuQjj ,  defined  in 
Eqs.  (3.28),  (3.31),  (3.39)  and  (3.40).  The  agreement  between  the  numerical  result 
and  our  statistical  model  for  the  different  cases  in  Fig.  3.6  shows  that  the  chaotic 
contribution  to  the  impedance  in  a  mixed  system  has  the  same  statistics  as  the 
impedance  in  a  purely  chaotic  system,  provided  one  accounts  for  variations  in  the 
size  of  the  chaotic  portion  of  phase  space  accessible  at  the  locations  of  the  ports. 

Our  second  goal  is  to  compare  the  previous  statistical  model  of  in  Ref.  [58, 
59]  (which  assumes  that  the  classical  trajectories  are  all  chaotic)  with  our  statistical 
model  of  Oj  [which  includes  chaotic  contributions  {£c,ij)  and  an  approximated  for¬ 
mula  for  regular  contributions  defined  in  Eq.  (3.42)].  Figure  3.7  shows  that 

our  statistical  model  (red  solid  curves)  predicts  the  probability  density  function  of 
much  better  than  the  previous  result  (blue  dashed  curves)  that  one  would  obtain 
by  supposing  that  the  entire  phase  space  was  chaotic. 

Note  that,  in  our  formulation  in  Eq.  (3.32),  4>mn(xi)  —  0  if  xt  is  located  in 
the  stem  of  the  mushroom.  Therefore,  if  at  least  one  port,  say  port  i,  is  located 
in  the  stem  of  the  mushroom,  then  =  0  and  only  chaotic  modes  contribute 
to  the  impedance,  i.e.,  ^  =  £c,ij-  In  the  insets  of  Fig.  3.7,  we  show  probability 
density  functions  of  ^R,tj  calculated  from  numerically  obtained  regular  eigenmodes 
and  the  probability  density  function  of  ^Rjij  calculated  from  our  approximate  regular 
eigenmodes  (delta  function  (red)  at  =  0  for  the  insets  to  Figs.  3.7(a  and  b)  and 
red  curve  in  the  inset  to  Fig.  3.7  (c)).  In  particular,  we  observe  that  the  pdf  widths 
in  the  insets  to  Figs.  3.7  (a  and  b)  are  much  less  than  for  the  inset  to  Fig.  3.7(c). 
The  small  pdf  widths  in  the  insets  to  Figs.  3.7  (a  and  b)  can  perhaps  be  explained 
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Figure  3.6:  Plot  of  the  probability  density  function  from  numerical  solution  (black 
histogram)  and  mean  probability  density  function  from  Monte  Carlo  simulation  (red 
solid  curve),  Eq.  (3.44),  with  root  mean  squared  error  bounds  (blue  dashed  curve), 
Eq.  (3.45).  The  black  and  red  dots  are  the  position  of  coaxial  transmission  lines 
(ports)  in  case  (a)  one  port  in  chaotic  region  and  the  other  in  mixed  region  (b)  both 
pots  in  chaotic  region  (c)  both  ports  in  mixed  region. 

by  dynamical  tunneling  (see  [80,  81]);  however,  this  effect  is  not  significant  in  the 
probability  density  function  of  ^  =  ^Rij  +  £c,ij  which  is  the  convolution  of  the 
probability  density  function  of  £c,ij  and  £,rtij- 
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(a)  (b) 


(c) 

Figure  3.7:  Plot  of  the  probability  density  function  of  ^  =  £7^  +  £c,ij  from  numer¬ 
ical  eigenmode  solution  (black  histogram),  our  statistical  model  that  treats  regular 
and  chaotic  contributions  separately  (red  solid  curve),  and  previous  statistical  model 
that  assumes  that  all  eigenmodes  are  chaotic  (blue  dashed  curve).  The  black  and 
red  dots  are  the  position  of  coaxial  transmission  lines  (ports)  in  case  (a)  one  port 
in  chaotic  region  and  the  other  in  mixed  region  (b)  both  pots  in  chaotic  region  (c) 
both  ports  in  mixed  region.  The  insets  show  the  probability  density  function  of  the 
regular  contribution,  for  numerical  eigenmode  solutions  (black  histogram)  and 
for  the  approximate  eigenmode  in  Eq.  (3.32)  (red  solid  curve). 
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3.4  Discussion 


In  this  chapter,  we  develop  a  method  for  obtaining  the  short  wavelength  sta¬ 
tistical  properties  of  the  impedance  matrix  of  wave  systems  whose  ray  equations 
yield  a  ‘mixed’  phase  space  with  coexisting  chaotic  and  regular  orbits. 

In  obtaining  our  results  for  the  mushroom  billiard,  we  assume  that  the  regular 
eigenmodes  are  approximately  the  same  as  the  eigenmodes  in  a  quarter  circle  cavity. 
In  formulating  our  theory,  we  have  neglected  the  possibility  that  there  may  be  some 
modes  where  the  regular  and  chaotic  phase  space  regions  are  coupled  by  dynamical 
tunneling,  thus  changing  both  the  eigenfunctions  and  eigenenergies.  These  mixed 
modes,  whose  eigenfunctions  show  characteristic  of  both  regular  and  chaotic  be¬ 
havior,  can  change  the  wave  scattering  properties  at  k2  near  these  resonances  and 
this  effect  can  be  treated  semiclassically  for  the  particular  modes  under  considera¬ 
tion.  However,  in  our  formulation,  we  are  not  interested  in  specific  k2  values  but 
rather  the  pdf  for  a  randomly  chosen  k2  values.  In  our  system  the  number  of  these 
chaos/regular  mixed  modes  appears  to  be  relatively  small  compared  with  modes  that 
are  predominantly  confined  to  either  the  regular  or  the  chaotic  phase  space  regions. 
Thus,  we  expect  mixed  chaos/regular  modes  do  not  make  a  significant  contribution 
to  the  mode  counting  formula  in  Eq.  (3.27),  and  this  expectation  is  confirmed  by 
the  good  agreement  between  our  numerical  results  and  theory. 

In  our  model,  appropriate  to  the  situation  that  we  numerically  tested,  we  as¬ 
sume  that  4>n(xi)  and  (j)n(xj)  are  independent  Gaussian  random  variables  for  chaotic 
wavefunctions,  which  only  applies  if  ports  i  and  j  are  far  apart,  k\xt  —  x3\  1,  and 
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both  ports  are  not  close  to  the  cavity  boundary.  This  assumption,  however,  is  not 
essential:  two-point  correlations  in  the  random  wave  model  have  been  previously 
studied  [78,  79]  and  can  be  accounted  for  by  regarding  4>n(xi)  and  4>n(xj)  as  corre¬ 
lated  bivariate  Gaussian  random  variables  with  a  correlation  that  takes  into  account 
direct  and  indirect  ray  paths  between  x%  and  x3  [82], 
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Chapter  4 


Conclusions  and  Future  Work 

In  chapter  ??,  we  have  presented  a  semi  classical  analysis  of  energy  level  split¬ 
ting  of  symmetric,  quantum-dot-type,  double-well  systems,  where  the  wells  are  sepa¬ 
rated  by  a  tunneling  barrier.  Our  analysis  quantitatively  explains  the  observed  mean 
splittings  and  their  fluctuations.  The  mean  is  found  to  be  independent  of  the  well 
shape  and  independent  of  whether  the  well  orbits  are  chaotic  or  not.  In  contrast, 
the  fluctuation  statistics  are  vastly  different  when  the  well  orbits  are  integrable,  as 
compared  to  when  they  are  chaotic,  with  the  chaotic  case  yielding  much  reduced 
fluctuations  when  the  lateral  barrier  length  is  large  compared  to  a  wavelength.  Fu¬ 
ture  elucidation  of  such  fluctuations  might  be  to  employ  Bogomolny’s  semiclassical 
Green  function  approach. 

In  chapter  3,  we  develop  a  method  for  obtaining  the  short  wavelength  statisti¬ 
cal  properties  of  the  lossless  impedance  matrix  of  wave  systems  whose  ray  equations 
yield  a  ‘mixed’  phase  space  with  coexisting  chaotic  and  regular  orbits.  We  use  a 
specific  mixed  system,  mushroom  billiard,  which  has  clearly  divided  phase  space. 
More  generic  systems  can  display  infinite  hierarchies  of  KAM  island  chains  encir¬ 
cling  other  KAM  island  chains  with  chaos  intermixed.  This  type  of  intermingling  of 
chaotic  and  nonchaotic  orbits  on  all  scales  is  not  present  in  the  mushroom  billiard 
where  non-smooth  shape  is  designed  to  yield  a  clear  division  between  chaotic  and 
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regular  phase  space  region.  Our  motivation  in  using  the  mushroom  shape  is  that  the 
simplicity  provided  by  its  clear  division  of  chaotic  and  regular  phase  space  allows  a 
potentially  simpler  theory.  We  hope  that  our  work  can  serve  as  a  basis  for  future 
study  applicable  to  the  case  of  generic  phase  space  structure. 
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Appendix  A 


Efficient  Way  to  Calculate  Eigenvalues  of  Gaussian  Random  Matrices 
In  each  realization  of  random  coupling  model  (RCM),  we  need  eigenvalues  of 
an  N  x  N  Gaussian  random  matrix,  where  N  >  103  in  general.  For  a  N  x  N  dense 
matrix,  general  eigensolvers  necessitate  0(N 3)  operations  to  solve  the  problem. 
Moreover,  we  also  need  more  than  1000  realization  to  get  a  good  ensemble  average,  so 
solving  the  eigenvalue  of  random  matrices  is  the  bottleneck  of  the  whole  calculation. 
In  this  appendix,  we  brief  the  tridiagonal  technique  that  allow  us  to  calculate  these 
eigenvalues  efficiently. 

Suppose  a  symmetric  N  x  N  Gaussian  orthogonal  random  matrix  H  =  (A  + 
At)/2,  where  A  is  also  an  N  x  N  matrix  with  AtJ  are  i.i.d.  N( 0, 1),  i.e., 

A(0,1)  A(0, 1/2) .  A(0, 1/2) 

N(  0,1/2)  •••  i 

H  ~  :  ••  ••  :  ,  (A.l) 

:  •••  N (0, 1/2) 

A(0, 1/2)  . A(0, 1/2)  A(0,1) 

J  NxN 

or 

{N( 0, 1)  for  i  =  j 

(A. 2) 

-A (0, 1/2)  for 

We  can  tridiagonalize  this  matrix  with  the  following  steps. 
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For  n  —  1,  we  define 


H«  =  H 


a'1’  =  -sgn(H^) 

\  3= 


N 


7=2 


r(1)  =  J ^(aW)2  -  H^laW 


v(i)  = 


n(1)lT 

7  UN  \ 


p^x7v  =  W-2v^(vW)^ 


(A.3) 

(A.4) 

(A.5) 

(A.6) 

(A.7) 


where 

( 


0 


< 


2  rW 


(l) 

fc,l 


2r(!) 


for  k  <  1 
for  k  =  2 
for  k  >  3  , 


(A.8) 


and  it  can  be  shown  with  some  algebra  manipulation  to  prove  that  P(1)  is 
symmetric  and  orthogonal,  and  thus  pd)Hd)p(1)  i§  an  orthogonal  transform 
and  we  can  define 


jj(2)  —  p(i)  jj(i)p(i) 


(A. 9) 


For  n  —  2, . . .. , N  —  2,  we  define 


N 

«(n)  =  sgn{H^lln) 

v 

(A-10) 

\ 

j=n+l 

r(n)  =  J^(aW)2-H, 

(">  a(") 

n+ljn^ 

(A. 11) 

v<">  =  [«<”>,  <>]T 

(A. 12) 

PnIn  =  I nxn  -  2v^(v^)T  (A. 13) 
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where 


( n )  . 

Vi—< 


H(hn) -q(") 

fc,n _ 

2r(n) 

(n) 


gj.,„ 

2r(n) 


for  k  <  n 
for  k  —  n  +  1 
for  k  >  n+2. 


(A.14) 


and  we  can  show  what  P(n)  is  symmetric  and  orthogonal;  thus,  P(n)H(n)p(n) 
is  an  orthogonal  transform  and  we  can  define 


jj(n+l)  _  p(n)jj(n)p(n) 


(A.  15) 


Finally,  we  reach 

JjOV-l)  _  p(AT— 2)  .  .  .  p(l)Hp(l)  .  .  .  p(Ar— 2) 


(A.  16) 


which  is  a  orthogonal  transform  of  H.  Thus,  the  eigenvalues  of  H  and  d 
should  have  same  distribution. 


Furthermore,  it  can  be  shown  that 


V2 


N(  0,2) 

Xn- i 

0 

0 

Xn- i 

N(  0,2) 

XN-2 

0 

Xn- 2 

0 

N(  0,2) 

Xi 

0 

0 

Xi 

N(  0,2) 

(A.  17) 


L  J  NxN 

where  Xn  is  a  Chi- distributed  random  variable  with  n  degree  of  freedom. 

For  a  symmetric  tridiagonal  NxN  sparse  matrix,  we  only  need  0(N2)  oper¬ 
ations  to  solve  the  eigenvalue  problem,  which  is  O(N)  faster  than  dense  matrix,  see 


Fig.  A.l. 
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N 


Figure  A.l:  Comparison  between  average  time  required  to  solve  the  eigenvalues  of 
a  N  x  N  GOE  matrix  using  (black)  Eq.(A.l)  and  (red)  Eq.(A.17)  in  log-log  plot. 
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Note  that,  this  method  apply  not  only  to  calculate  the  eigenvalues  of  GOE 


matrices  but  also  GUE  and  GSE  matrices. 


H 


(iV-l) 


1 

71 


^(0,2)  X/3(n- i)  0  ...  0 

X/3{N-1)  N(0,2)  X/3(N-2)  '•  : 

o  X/3(N- 2)  '  •  '  •  0 

:  -V(().  2)  Xp 

o  ...  o  xp  #(0,2) 


NxN 


where  [3—1  for  GOE  matrices,  [3  —  2  for  GOE  matrices  and  (3  = 
matrices. 

Matlab  source  code  for  Eq.  (A.  18) 

°/c  off-diagonal  term 

d0  =  normrnd(0 , sqrt (2)  , 1  ,N)  1  ; 

°/,  diagonal  term 

dl=sqrt(chi2rnd(beta*((N-l) : - 1 : 1 ) ) ) ’ ; 

°/,  form  tri-diagonal  matrix 

H=spdiags (  [  [dl ; 0]  ,d0,  [  0 ; d 1 ] ]  ,  [-1,0,1]  ,N,N) /sqrt (2) ; 

%  solve  eigenvalue 
ks=eig (H) ; 


(A.18) 


for  GSE 
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Appendix  B 


Method  of  Particular  Solutions 

In  chapter  3,  numerical  calculation  of  eigenvalue  problems  of  the  two-dimensional 
Helmholtz  equation  plays  the  major  role  to  examine  our  theoretical  prediction.  In 
this  appendix,  we  briefly  introduce  numerical  techniques  to  solve  the  two-dimensional 
Helmholtz  equation  efficiently  and  more  detail  can  be  found  in  [63,  73,  71,  83].  Nu¬ 
merical  result  of  two  billiards  (stadium  and  mushroom)  has  been  provided.  New 
discover  of  over  counting  eigenmodes  in  previous  algorithm  has  been  discussed. 
More  statistical  properties  of  eigenvalues  and  eigenfunctions  will  be  covered  in  Ap¬ 
pendix  C. 


B.l  Introduction 


Considering  w(x)  satisfies  the  Helmholtz  equation  with  Dirichlet  boundary 
condition  in  a  two-dimensional  domain  B 


(V2  +  k2)u(x )  =0  for  x  G  B 
■u(x)  =0  for  x  G  <9B, 


(B.l) 

(B-2) 


with  normalization  condition 


(B.3) 


83 


B.2  Method  of  Particular  Solutions 


The  idea  of  method  of  particular  solutions  (MPS)  is  to  approximate  the  eigen¬ 
function  with  wavenumber  k ,  u(k]x),  by  a  linear  combination  of  basis  function, 
{0i (k;  x),4>2 (k]x), . . which  satisfy  Eq.  (B.l)  but  not  necessary  satisfy  the  bound¬ 
ary  condition  (B.2),  i.e., 

u(k-x)  =  y]cm0m(fc;£).  (B.4) 

m 

We  dehne  boundary  tension  function  and  area  norm  for  u(k\  x)as 


Tb[u]  =  (j  \u(k;  s)\2  ds, 

(B.5) 

Jan 

Tj[u]  =  (  \u(k;  x)  2d2x. 
in 

(B.6) 

If  kn  is  an  eigenvalue  of  Eqs.  (B.l)  and  (B.2),  we  can  find  a  nonzero  set  of  coef¬ 
ficient,  {cnm},  such  that  TB[u }  =  0.  Moreover,  u(kn;x )  should  always  satisfies  the 
normalization  condition,  i.e.,  Tj[u]  =  1. 

To  solve  kn  and  {cnm},  we  discretize  Eqs.  (B.5)  and  (B.6)  as 

Mb 

Tb[u\  =  ^2  \u(xb)\2Asb,  (B.7) 

b=lixj)Gdfl 

Mj 

Ti[u\  =  \u(xi)\2  Avi.  (B.8) 

i=l,Xi^Q 

For  simplicity,  we  choose  A sb  =  As  for  all  sb  and  A <ji  =  Aa  for  all  at.  Finding 
zero  of  T B [u]  usually  reaches  a  trivial  solution,  cmn  =  0,  for  all  m.  To  avoid  this 
undesired  solution,  we  try  to  find  zeros  of  Tb[u]/Tj[u\;  however,  we  still  cannot  solve 
{him})  since  there  are  infinity  terms  of  {cnm},  so  we  truncate  the  infinite  sum  in 
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Eq.  (B.4)  by  a  finite  sum 


M 

u(kn;x)  f»  X  C-nm^m,  (kn;x).  (B.9) 

m 

Instead  of  finding  zeros  of  Tb[u]/Tj[u],  we  look  for  local  minimums  of  TB[u]/Tr[u], 
We  define  a  MB  x  M  matrix,  AB{k),  and  a  Mj  x  M  matrix,  A i(k),  and  their 
matrix  elements  are 


AB,bm  =  4>m(xb),  for  xh  £  <90,  6=1,...,  Mb  (B.10) 

for  £  0,  i  =  1, . . . ,  Mj  (B.ll) 

Thus,  finding  local  minimums  of  T#  [w] /Tj  [w]  is  equivalent  to  finding  local  minimum 
of 


(B.12) 


TY[u]  lA/^cl2’ 

where  c  =  [ci, . . . ,  cm]t,  |  •  |  is  the  Euclidean  norm,  and  this  problem  can  be  solved 
by  the  general  singular  value  decomposition  (GSVD). 

Assuming  A  e  9ftnxp  and  B  e  3fJmxp  and  n  >  p.  There  exist  orthogonal 
matrices  U  £  Krixn  and  W  £  9£mxm  and  an  invertible  matrix  X  £  such  that 


A  =  UCX1,  B  =  WSX  \ 

where  C  £  3fJnxp  and  S  £  ffinxP  are  diagonal  matrices  with  0  <  C\  <  ■  ■  ■  <  cp  <  1  and 
0  <  Sj  <  1  for  j  —  1, . . . ,  min {m,p}  with  sj  +  c|  =  1  and  Cj  =  1  for  j  >  min {m,p}, 
i.e.,  Sj  =  0  for  j  >  mm{rn,p} .  The  values  a3  =  Cj/sj  are  the  generalized  singular 
values  of  A  and  B.  Define  x*  =  [Au, . . . ,  Xpj]T,  we  have 

|  Ax?:  | 2  =  c2 

|  Bx,;  | 2  =  s2.  (B.13) 
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Applying  the  GSVD  to  our  problem,  we  assume  that  A  =  A B(k)  and  B  = 
A i(k)  and  choose  number  of  boundary  discrete  point  greater  than  basis.  We  do  not 
describe  the  procedure  of  GSVD  here,  but  in  MATLAB,  there  is  a  built-in  gsvd() 
function  that  allow  us  to  perform  the  GSVD  easily.  At  each  k,  we  can  find  out  the 
minimal  singular  value,  say  (J\{k)\  scanning  over  k,  we  can  find  out  at  kn ,  o\ (kn)  is 
the  local  minimal.  From  Eq.  (B.13),  we  know  the  corresponding  cn  =  y^i{kn). 

The  method  of  particular  solutions  can  solve  the  eigenvalue  problem  precisely 
by  scanning  over  k  very  carefully,  but  slow  down  the  performance;  furthermore, 
this  method  easily  neglects  one  of  two  eigenmodes  whose  eigenvalues  are  close  each 
other.  In  practical,  we  only  use  this  method  to  find  out  the  first  few  eigenmodes. 


B.3  Scaling  Method  for  MPS 


Given  a  wavefunction  ^(/yr),  we  can  define  its  scaling  function  near  a  given 
wave  number,  /c0,  as 

ip(k]r)  =  -0  r  (B.14) 

the  first  derivative  over  k  near  ko  is 


a=l 


=  —  r  ■  VMar), 
k0 


(B.15) 


and  the  Taylor  expansion  for  the  wavefunction  can  be  written  down  as 


P(^o  +  8',r)  = 


1  +  V)  +  \m{?'  V)2 + 0(,53) 


(B.16) 


For  s  E  dkl,  since  V  =  hdn  +  tdt,  where  t  and  n  are  tangential  and  outgoing  normal 
direction  at  boundary,  see  Fig.  B.3(a),  we  can  write  down  the  scaling  function  for 


an  exact  eigenfunction,  ip^k^r),  with  (V2  +  fc2)V;/j  =  0  and  Dirichlet  boundary 
condition,  at  k  =  +  8  as 


\  6  8 2  82  a 

+  8;  S)  =  -j-rnd„  +  -T^rnrtdndt  +  T^--(r2  -  r2n)dn 

where  a  is  the  inverse  of  radius  of  curvature  at  s. 


^(s)  +  0(83),  (B.17) 


For  a  general  sum  of  scaling  eigenfunction,  i.e., 


^(/c;r)  =  (B.18) 

the  tension  is 

f{k)=  (l  dsw(s)\i/}(k]s)  |2.  (B.19) 

JdD 

where  w(s)  =  l/rn  is  a  boundary  weighting  function.  The  tension  can  be  expressed 
as  a  quadratic  form 

f(k)  =  xTF(A;)i,  (B.20) 


where  the  matrix  element  of  F(k)  is 


Ffjvik)  =  (b  dsw(s),ijjfl(k1  s)ipv(k,  s) 
JdD 

Sv  8V 


kfj,  kv  j  Q£) 


=  TT~  f  s))(dn^u(k,  s))  +  0(83) 


=  2  S^M^ik)  +  0(83) 


(B.21) 


and 


M^ik)  =  —  j)  dsw^r^dn^^k,  s))(dni>u(k,  s)) 


(B.22) 


where  we  assume  that  /qt  —  ku  =  k. 


It  has  been  proven  that  [83] 


Mfj,u(k)  fa  8^. 
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(B.23) 


Combining  with  Eq.  (B.21)  and  using  the  fact  that  =  1  [since  k  =  +  5],  we 

get 

dF 

~^-{k)  =  2  {5,  +  SU)M,U  +  0(52).  (B.24) 

Now,  we  want  to  express  the  scaling  eigenfunction  as  a  linear  combination  of 
basis  function,  i.e. , 

M^)  =  (R25) 

i 

and  we  have 

F  (k)  =  YTF(k)Y,  (B.26) 

fW=Y-fWY,  (B.27) 

where  the  matrix  element  of  F (k)  and  <9F (k)/dk  are 

Fij(k)  —  ®  ds — 4>i(k;  s)(ftj(k]  s)  (B.28) 

JdD  r  n 

—Tr~(k)  =  T(£  ds — (pAk;  s)s  ■  Y4>j(k;  s)  +  transpose.  (B.29) 

dk  k  JQD  r n 

Combining  Eqs.  (B.21),  (B.24),  (B.23),  (B.26),  (B.27),  (B.28)  and  (B.29),  we  can 
solve  the  following  generalized  eigenvalue  problem 

(d-A'‘F)y'‘  =  °  <R3°) 

where  =  j-  is  the  generalized  eigenvalues  such  that  k^  =  k  —  5^  and  = 
[Yim,  . . . ,  YMfJ\T  is  the  generalized  eigenvectors  such  that  ^(r)  = 

The  advantage  of  scaling  method  for  MPS  is  that  at  each  time  we  solve  the 
generalized  eigenvalue  problem  of  F  and  dF / dk ,  Eq.  (B.30),  we  would  be  able  to  solve 


all  eigenmodes  near  k.  Using  proper  number  of  basis  function  and  discrete  boundary 


point  (see  Refs[  add  something]),  we  would  be  find  out  all  eigenmodes,  E  [k  — 
dk,  k+dk]  by  solving  the  generalized  eigenvalue  problem  once.  To  solve  all  eigenvalue 
k/t  <  kMax,  we  only  need  to  choose  proper  dk  and  solve  N  generalized  eigenvalue 
problems  on  the  following  N  interval,  [0,  2dk],  (2dk,  Adk], . . . ,  (2(N  —  l)dk,  2Ndk], 
where  2 (N  —  1  )dk  <  kMax  <  2Ndk. 

B.4  Over  Counting  Eigenmodes 

The  scaling  method  for  MPS  allows  us  to  find  out  all  eigenmodes  of  two- 
dimensional  Helmholtz  equation  with  Dirichlet  boundary  condition.  Actually,  it 
find  out  more  eigenmodes  than  the  Weyl’s  formula  prediction.  To  explain  this  over 
counting  phenomena,  we  focus  on  all  eigenmodes  we  solved  in  the  two  consecutive 
interval,  ( k  —  dk,  k  +  dk]  and  ( k  +  dk,  k  +  3  dk],  say  k  —  dk  <  ...  <  ...  <  kn_\  <  kn  < 
k  +  dk  and  k  +  dk  <  kn+ 1  <  kn+ 2  <  . . .  <  k  +  3dk.  Since  the  numerical  solution  is  not 
absolutely  precise  to  the  exact  eigenvalue,  it  is  possible  that  kn  and  kn+ 1  actually 
represent  the  same  eigenmodes.  To  eliminate  these  over  counting  eigenmodes,  we 
can  solve  the  generalized  eigenvalue  problem  in  ( k ,  k  +  2 dk]  to  verify  whether  there 
are  one  or  two  eigenmodes  near  k  +  dk  (see  Fig.  B.l). 

We  compare  the  difference  mode  counting  function  from  numerically  calculated 
eigenmodes  of  stadium  billiard,  see  Fig.  B.3(a),  and  the  Weyl’s  formula  (1.8)  in 
Fig.  B.2,  and  mode  counting  function  excluding  the  over  counting  eigenmodes  agrees 
with  Weyl’s  formula  better. 
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kn—1 


kn+1 

kn+2 


Figure  B.l:  Illustration  of  over  counting  eigenmodes. 


Figure  B.2:  Differencee  between  mode  counting  function  from  numerical  calculated 
eigenmodes  (black  solid  curve)  including  over  counting  eigenmodes  (red  solid  curve) 
excluding  over  counting  eigenmodes  and  Weyl’s  formula  prediction.  The  vertical 
blue  dashed  line  label  the  location  of  over  counting  eigenmodes. 

B.5  Choice  of  Basis 

The  most  natural  choice  of  basis  is  motivated  by  the  random  plane  wave 
hypothesis.  We  use  the  stadium  billiard  as  our  example,  see  Fig.  B.3(a).  We  can 
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choose  r0  G  fl,  and 


cos  [km  •  (r  —  r0)]  for  n  =  2 m  —  1 

sin  [km  •  (r  —  r0)]  for  n  =  2m 


Mr)  = 


(B.31) 


where  km  =  k  cos  (^)  x  +  k  sin  (^)  y,  n  —  1, . . . ,  2 M.  Some  eigenfunctions  are 
plotted  in  Figs.  B.4  and  B.5  . 

In  case  that  has  a  corner,  like  the  mushroom  billiard(  B.3(b)),  we  can  choose 
Fourier  Bessel  basis  function 


Mr)  =  Jan{kr )  sin  (and),  (B.32) 

where  Jan  is  the  Bessel  function  of  the  first  kind  with  order  an.  In  our  numerical 
work,  we  choose  a  =  3/4.  Some  eigenfunctions  are  plotted  in  Figs.  B.6  and  B.7  . 


Figure  B.3:  (a)  Stadium  billiard  and  (b)  Mushroom  billiard. 
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Appendix  C 


Statistics  of  Eigenmode 

In  this  section,  we  compare  the  numerically  calculated  wavefunction  amplitude 
with  the  prediction  of  random  matrix  theory  and  random  plane  wave  hypothesis. 
We  use  the  numerically  calculated  eigenmodes  of  stadium  and  mushroom  billiard  in 
Appendix  B  as  our  example  of  chaotic  and  mixed  system. 

C.l  Chaotic  Billiard 

We  solve  the  first  23,072  eigenmodes  (kn  <  300)  of  the  stadium  billiard,  see 
Fig.  B.3(a),  use  Eq.  (1.11)  to  normalize  the  spacing  between  two  consecutive  eigen¬ 
value  and  get  good  agreement  with  the  random  matrix  theory  (1.14),  see  Fig.  C.l. 

We  also  compared  the  eigenfunction  amplitude  at  different  position,  0(x),  in 
the  billiard  and  compare  the  probability  density  function  of  0(x)  with  Eq.  (1.18) 
and  get  good  agreement,  see  Fig.  C.2. 
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Figure  C.l:  Comparison  between  the  probability  density  function  of  normalized 
nearest  neighbor  eigenvalue  spacings  of  (black  histogram)  stadium  billiard  and  (red 
curve)  Wigner  GOE  distribution. 

C.2  Mixed  Billiard 


We  solve  the  first  39,114  eigenmodes  ( kn  <  500)  of  the  mushroom  billiard,  see 
Fig.  B.3(b),  and  use  the  outward  normal  derivative  of  eigenfunction  at  the  bound¬ 
ary  to  classify  whether  eigenmodes  are  either  regular  (0r,  k2)  or  chaotic  (0C,  k2), 
see  Sec.  3.3.  In  Fig.  C.3,  we  use  Eq.  (1.11)  to  normalize  the  spacing  between  two 
consecutive  ( chaotic/regular /mixed)  eigenvalue  and  get  good  agreement  with  the 
random  matrix  theory  (1.14),  Poisson  distribution  (1.16),  and  Berry-Robnik  distri¬ 
bution  (C.l), 

d2 

Pbr{Pv,  Pc,  s)  — 

where 

2  f°° 

erfc(x)  =  —  exp  (—t2)dt.  (C.2) 

V  Jx 


exp  (— prs)erfc  f  pcs 


7 T 


(C.l) 
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Figure  C.2:  Comparison  between  (black  histogram)  the  probability  density  func¬ 
tion  of  the  first  23,072  eigenfunction  amplitudes  of  the  stadium  billiard  at  different 
position  (see  inset)  and  (red  curve)  Gaussian  distribution. 


Note  that  in  Berry- Robnik  distribution  (C.l),  it  requires  the  information  of  phase 
space  volume  ratio  of  regular  and  chaotic  region,  i.e.,  pr  and  pc.  Using  Eqs.  (3.28) 
and  (3.37),  we  get  pc  =  AJA  and  pr  =  1  —  pc,  where  A  is  the  total  area  of  the 
mushroom  cavity. 

We  also  compared  the  chaotic  eigenfunction  amplitude  at  different  position, 
</>c(x)>  in  the  billiard  and  compare  the  probability  density  function  of  0c (x)  with 
Eqs.  (3.30)  and  (3.31)  and  get  good  agreement,  see  Fig.  C.4. 

There  is  no  universal  theorem  to  predict  the  probability  density  function  of 
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Figure  C.3:  (a)  Comparison  between  the  probability  density  function  of  normalized 
nearest  neighbor  chaotic  eigenvalue  spacings  of  (black  histogram)  mushroom  billiard 
and  (red  curve)  GOE  matrix,  (b)  Comparison  between  (black  histogram)  the  prob¬ 
ability  density  function  of  normalized  nearest  neighbor  regular  eigenvalue  spacings 
of  mushroom  billiard  and  (red  curve)  Poisson  distribution,  (a)  Comparison  between 
(black  histogram)  the  probability  density  function  of  normalized  nearest  neighbor 
eigenvalue  spacings  of  mushroom  billiard  and  (red  curve)  Berry-Robnik  distribution. 

0#(x)  over  a  wide  range  of  k 2,  see  Fig.  C.5. 
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Figure  C.4:  Comparison  between  (black  histogram)  the  probability  density  function 
of  the  first  19,198  chaotic  eigenfunction  amplitudes  of  the  mushroom  billiard  at 
different  position  (see  inset)  and  (red  curve)  Gaussian  distribution. 
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Figure  C.5:  The  probability  density  function  of  the  first  19,916  regular  eigenfunction 
amplitudes  of  the  mushroom  billiard  at  different  position  (see  inset). 
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Appendix  D 


Random  Wave  Model  with  Boundary  Conditions 

Given  a  two-dimensional  chaotic  region  hi  with  area  A,  we  assume  that  ri;  e 
hi  and  |  r7;  —  r^j  much  less  than  the  distance  of  r*  or  r j  to  the  boundary  dfl. 
From  Berry’s  conjecture,  the  eigenfunction  amplitudes  at  r*  and  rj,  (0(rd,  <P(*j))  = 
(c f>i,(j)j ),  are  normal  distributed,  and  the  correlation  between  A,  A  is  Jo(k\r-ir _,j). 
Thus,  ((pi,  <pj)  can  be  described  as  a  bivariate  Gaussian  distribution 


Gsd(ri?rp  k)  =  - 


ih(2'nih)^d~1^2 


£ia, 


1/2  ^iSp/h—ii/piT  /  2 


(D-3) 


where  the  sum  is  over  all  classical  paths  connecting  between  and  r,-,  vp  is  number 


of  classical  focal  points  along  the  path, 


(D.4) 
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is  the  actions  of  the  path  and 


D 


p 


d2nSp 
dr  j dr , 


d2nSp  \ 
dEdr.i  I 


(D.5) 


I  d2nSp  d2nSp  I 
y  dr  jE  dE'2  J 

is  the  determinant  of  the  second  derivative  matrix  of  Sp.  Then  the  two-point  corre¬ 
lation  function  at  a  given  wavenumber  k  is 


C(ThTj\k)  = 


2n  ip(k2 


[Gsci(rj,  r p  k)*  -  Gsd(rj,  rp  k)] 


p(fc2)(2^,)(rf+1)/2  S  W2  cos  ~  (2ur  +  d~  1)7r/4]'  (D-6) 


paths 

Using  the  asymptotic  form  of  the  Bessel  function  of  the  first  kind, 


Mz)  “  V^cos  (z  “  i)  • 


(D.7) 


we  get 

1 

C0(ri,  r j,  k )  =  —J0(k\Ti  -  | ) , 

which  is  identical  with  Berry’s  original  result. 

Taking  r,:  =  rv  we  get 


cr2(r  p  k )  =  C(rj,  rj,  k)  =  —  [1  T  boundary  term  (A;)] , 

yjL 


(D.8) 


(D.9) 


which  is  depend  on  wavenumber  k  and  position  r*.  To  examine  this  boundary  effect, 
we  average  over  k 2  from  0  to  k2M 


/  27  J*M  C(ri,ri,k)p(k2)dk2 

(a  (rj)  =  -u 


foM  p(k2)dk2 


(D.10) 


and  replace  <r2(x)  in  Eq.  (3.30).  In  Fig.  D.l,  the  semiclassical  correction  predict 
the  average  eigenfunction  density  over  k2  agree  with  numerical  calculated  23,  072 
stadium  eigenfunctions  pretty  well 
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(a) 


r 

(b) 

Figure  D.l:  (a)  Stadium  Billiard  (b)  Comparison  between  (Black  Dots)  average  over 
n  —  1, . . . ,  23072  of  numerical  calculated  <^(r)  and  (red  curves)  Eq.  (D.10). 
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Appendix  E 


Lorentzian  Distribution  of  the  Regular  Normalized  Impedance 

Consider  the  normalized  impedance, 


N 


-*3 


1  \  WniWnj 

n  ^[~k2  ~k2' 


(E.l) 


where  (wni,  wnj)  are  bivariate  random  variables  with  probability  density  function 
(PDF)  fij(wni,wnj),  and  k2  are  independent  random  variables  distributed  uniformly 
on  (0 ,k%),  i.e.,  the  PDF  is  /^(k2)  =  1  /k2N.  Let 


*  _  1  U-'ni^ri'j 

7  Q  T  O  ^ 

vr  k2  -  k2 


(E-2) 


such  that 


“ ij  —  £W  (E-3) 

n 

The  PDF,  /h(^),  and  the  characteristic  function  of  5,  $=(£)  are  given  by 

N 

Hz)  =  /  (e.4) 

^  n= 1  n' 

N 

$=(7)  /  d^i . .  .d^M  |  A  (A)  exp  (if  E«»0  =  [*5(*)]W.  (E-5) 

n=  1  n! 

where  A(A)  is  the  PDF  of  £  and  ^(f)  =  J  d£n  exp(it£n) /t(£n)  is  the  characteristic 
function  of  £.  We  can  calculate  $^(t)  by  directly  evaluating  the  integral, 


<f>?(t)  =  /  dwnidwnjfij(wni,w, 


nj) 


[kN  _  1 

x  /  dki  — —  exp 

Jo  nk% 


■  1 1  UZniWrij 

1  ^k2-k2 


(E.6) 
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For  small  values  of  t,  relevant  in  the  limit  N  3>  1,  the  second  integral  of  (E.6)  is 


[Nd~enew 

k2N  Jo  \  tt  k2  -  k2 


=1  + 


U'njVJnj  .  1  U-niWyij 


_  _  +0(t2), 

jr  U‘2  p  —  h2  ^ 


(E-7) 


which  to  first  order  in  t  yields 


$t(t)  ~  /  dwnidwnjfij(wni,wr 


(\t\\wniWnA  1  WniWni  ,  k2 

1  +  1  "  -  -  ft - log  ^ 

U2  7f  U2  b.2  _  L.2 

rvjy  rvjy  tv 

_ ;j_E{wnjWnj}  k 2 

k\  V  71  k%-k2 


I  ^  |  E  { |  wniWnj  i} . 


where  E{- ■  =  /  •  •  •  fij(wni,  wnj)dwnidwnj.  Now,  we  calculate  $=(t);  since  the 

mean  spacing  between  adjacent  k2  is  normalized  to  unity,  we  can  replace  k2N  in  (E.8) 
by  N  and  insert  it  into  (E.5).  As  N  — y  oo,  we  obtain 

~  U\  [i  1  (  -  E  {WniWnj}  k2 

k«i^ 


1 1\  E  { \  wniWnj  | } 


(  E{wniwnj}  k 2 

->•  exp  it - log  -  \t\E{\wniwnj\} 

\  ^  kM  fc 


Comparing  with  the  characteristic  function  of  a  Lorentzian  RV  with  mode  xq  and 
width  W,  <h(t)  =  exp  (itx0  —  W|t|),  we  know  is  Lorentzian  distributed  with 
mode  E {wniWnj }( log  \k2\  —log  | k%  —  k2\)/7t  and  width  E{\wniWnj\}.  Since  the  spac¬ 
ing  distribution  of  k2  for  regular  systems  is  exponential  distributed,  as  N  — >  oo,  the 


distribution  of  k2  is  uniformly  distributed  in  (0,  N );  thus,  the  normalized  impedance 
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of  regular  systems  are  also  Lorentzian  distributed  and  all  the  system  specific  infor¬ 
mations  are  included  in  mode  and  width  of  the  Lorentzian. 
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